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TURBULENCE  IN  COMPRESSIBLE  FLOWS 

Final  Report  for  Grant  F49620-97- 1-0089 

Robert  D.  Moser 

Department  of  Theoretical  and  Applied  Mechanics 
University  of  Illinois,  Urbana-Champaign 


1  Motivation  and  Objectives 


A  large  fraction  of  the  effort  in  turbulence  research,  including  numerical  simulations,  ex¬ 
periments,  theory  and  modeling,  has  been  directed  toward  incompressible  turbulence.  This 
is  largely  due  to  the  simplifications  that  the  incompressible  assumption  allows,  and  to  the 
belief  that  turbulence  is  in  essence  an  incompressible  phenomenon.  However,  many  of  the 
flows  of  technological  interest  are  compressible  (e.g.  external  aerodynamics,  propulsion 
systems)  leading  to  the  requirement  that  our  understanding  of  incompressible  turbulence 
be  translated  to  these  compressible  flows.  An  example  of  this  is  that  most  commonly  used 
turbulence  models  for  compressible  flows  are  modifications  of  models  developed  for  in¬ 
compressible  flows  (e.g.  the  k-e  model). 

Fortunately,  in  many  flows  the  small-scale  turbulence  is  nearly  incompressible,  albeit  with 
non-constant  density,  even  though  the  flow  itself  is  highly  compressible.  For  example,  in 
much  of  the  literature  on  compressible  boundary  layers  (see  the  review  by  Spina  et  al., 
1994),  differences  relative  to  an  incompressible  boundary  layer  (up  to  M  =  5,  say)  are 
understood  as  being  caused  by  variations  in  the  mean  density  across  the  layer.  This  has  long 
been  appreciated  and  is  the  basis  of  the  Van  Driest  (Van  Driest  1951, 1956)  transformation 
and  the  Morkovin  hypothesis  (Morkovin  1961). 

These  observations  suggest  that  a  formulation  based  on  low  Mach  number  asymptotics,  in 
which  the  turbulence  is  treated  as  weakly  compressible  while  the  mean  is  considered  to  be 
fully  compressible,  would  provide  a  good  description  of  the  boundary  layer  to  quite  large 
Mach  numbers,  and  a  rational  basis  for  the  development  of  compressible  turbulence  mod¬ 
els.  Aided  by  direct  numerical  simulation  data,  the  research  under  this  grant  is  aimed  at 
developing  and  evaluating  such  weakly  compressible  descriptions  of  compressible  turbu¬ 
lence. 


2  Supported  Research 


There  were  four  primary  research  activities  supported  under  this  grant,  all  aimed  at  the 
weakly  compressible  analysis  of  turbulence  in  compressible  shear  layers.  These  activities 
are  listed  below: 


1.  Direct  numerical  simulation  of  a  compressible  boundary  layer  at  Ma  =  2.5 

2.  Low  Mach-number  asymptotic  analysis  and  numerical  decomposition  of  compress¬ 
ible  boundary  layers. 

3.  Development  and  implementation  of  a  new  B-spline  collocation  numerical  method 
for  the  numerical  simulation  of  compressible  shear  layers. 

4.  Numerical  simulation  and  analysis  of  compressible  mixing  layers. 

The  key  results  of  each  of  these  activities  will  be  described  briefly  below.  More  details  are 
available  in  the  papers  attached  as  appendices. 


2.1  DNS  of  a  Compressible  Boundary  Layer 

Using  a  technique  similar  to  that  of  Spalart  (1989)  to  homogenize  the  streamwise  direc¬ 
tion,  a  Mach  2.5  compressible  boundary  layer  was  simulated.  A  B-spline  Galerkin  method 
was  used  and  the  simulation  was  performed  in  a  computational  domain  large  enough  to 
eliminate  most  finite  domain  size  effects.  The  simulation  was  used  to  evaluate  a  variety  of 
approximations  used  in  the  modeling  of  compressible  boundary  layers.  For  example,  the 
Strong  Reynolds  Analogy  (Morkovin,  1961)  was  evaluated  and  found  to  be  in  disagree¬ 
ment  with  the  DNS  data  when  applied  in  its  most  common  form.  Valid  approximations 
were  also  determined.  Of  particular  interest  to  this  project  was  the  determination  that  the 
velocity-temperature  correlation  was  consistent  with  that  in  an  incompressible  boundary 
layer. 

The  details  of  the  simulation  and  the  results  were  reported  in  Guarini  et  al  (2000),  a  reprint 
of  which  is  included  as  appendix  A.  This  research  activity  was  also  sponsored  by  the 
NAS  A- Ames  Research  Center  through  a  Graduate  Research  Fellowship  for  Guarini.  The 
efforts  of  the  PI  and  computer  resources  at  HPCMP  were  supported  through  the  current 
grant. 


2.2  Analysis  of  Compressible  Turbulent  Boundary  Layers 


A  weakly  compressible  asymptotic  analysis  similar  to  that  of  Zank  &  Matthaeus  (1991) 
was  developed  for  the  turbulence  in  a  compressible  boundary  layer.  In  developing  such 


an  asymptotic  analysis,  a  variety  of  different  assumptions  can  be  made  regarding  scaling 
of  different  components  of  the  compressible  fields.  The  direct  numerical  simulation  data 
of  Maeder,  Adams  &  Kleiser  (2000)  at  three  Mach  numbers  (3,  4.5  and  6)  were  analyzed 
to  determine  the  appropriate  scalings.  Further,  it  was  determined  that  unlike  the  heat-flux 
dominated  approximations  of  Zank  &  Matthaeus  (1991),  the  acoustic  contributions  to  some 
quantities  (e.g.  pressure  fluctuations  and  velocity  dilatation)  are  not  negligible.  To  evaluate 
the  impact  of  acoustics  on  the  turbulent  boundary  layer,  a  decomposition  procedure  was 
devised  to  allow  a  compressible  turbulent  field  to  be  decomposed  into  acoustic  and  nona¬ 
coustic  parts.  This  decomposition  was  applied  to  the  DNS  data  of  Maeder  et  al  (2000).  The 
asymptotic  analysis  suggests  and  the  decomposition  analysis  confirms  that  the  acoustic 
fluctuations  are  only  weakly  coupled  to  the  nonacoustic,  even  at  Ma  =  6.  The  nonacoustic 
equations  are  thus  autonomous,  and  can  be  considered  to  be  the  governing  equations  for 
the  turbulence  in  the  compressible  boundary  layer.  These  equations  can  thus  be  used  as  a 
basis  for  further  simulation  and  modeling. 

The  decomposition  analysis  and  its  results  are  described  in  a  paper  under  consideration  for 
publication  in  Theoretical  and  Computational  Fluid  Dynamics  (Borodai  &  Moser,  2000), 
and  a  preprint  is  included  as  Appendix  B. 

The  consequences  to  modeling  approximations  of  the  low  Mach  number  asymptotics  and 
the  autonomous  nonacoustic  governing  equations  is  being  explored.  Several  are  already 
clear.  For  example,  this  analysis  shows  that  the  pressure-dilatation  correlation  and  the 
so  called  dilatational  dissipation,  both  terms  in  the  equation  for  turbulent  kinetic  energy, 
should  be  negligibly  small,  as  has  been  recently  observed  in  DNS  data.  The  analysis  also 
suggests  that  the  statistics  of  the  compressible  turbulent  boundary  layer  should  be  related  to 
those  of  incompressible  boundary  layers  with  passive  scalars.  This  connection  is  continuing 
to  be  explored. 


2.3  Development  of  Compressible  Shear  Layer  Code 

Experience  with  the  boundary  layer  simulation  described  above  indicated  that  the  Galerkin 
formulation  had  several  practical  shortcomings,  the  most  obvious  being  inefficiencies  in¬ 
herent  in  evaluation  of  nonlinear  terms.  B-spline  collocation  methods  do  not  suffer  from 
this  shortcoming.  Furthermore  it  was  determined  that  they  have  several  other  desirable 
features.  In  particular,  these  collocation  methods  have  better  resolution  properties  than  the 
related  Galerkin  methods.  A  detailed  comparison  of  these  B-spline  collocation  methods 
with  other  high-resolution  methods  was  conducted,  and  as  a  result,  they  were  selected  as 
the  basis  for  a  new  implementation  of  a  compressible  shear  flow  code.  This  code  is  being 
used  for  the  mixing  layer  simulation  discussed  in  2.4  below. 

The  results  of  the  numerical  methods  evaluations  are  the  subject  of  a  paper  being  considered 
for  publication  in  Journal  of  Computational  Physics  (Kwok,  Moser  &  Jimenez,  2000),  and 
a  preprint  of  this  paper  is  included  as  Appendix  C. 


Figure  1:  Rms  of  a)  velocity  divergence  and  b)  temperature  fluctuations  in  a  two- 

dimensional  Mac  —  0.4  mixing  layer.  Shown  are - total, - ,  nonacoustic  and 

- acoustic  contributions 

2.4  Simulation  and  Analysis  of  Mixing  Layers 


The  compressible  shear  layer  code  described  in  2.3  is  being  used  to  perform  simulations 
of  compressible  mixing  layers.  Free  shear  layers,  such  as  the  mixing  layer  are  known  to 
exhibit  much  stronger  compressibility  effects  than  boundary  layers  at  similar  free-stream 
Mach  numbers.  They  are  thus  of  great  interest  in  investigating  the  coupling  of  acoustic 
fluctuations  with  turbulence.  A  decomposition  analysis  similar  to  that  for  the  boundary 
layer  is  being  applied  to  the  mixing  layer. 

Two-dimensional  mixing  layer  simulations  with  convective  Mach  number  of  0.4  and  0.8 
were  carried  out  to  provide  input  for  both  the  theory  development  and  planned  three  di¬ 
mensional  simulations.  In  particular,  the  simulations  provide  critical  information  on  the 
relative  importance  of  incompressible,  acoustic  and  thermal  parts  of  the  velocity,  tempera¬ 
ture,  pressure  and  density  fields.  For  a  mixing  layer  with  equal  free-stream  temperatures, 
we  expect  the  thermal  fluctuations  to  be  less  important  than  in  the  boundary  layer,  since 
the  mean  temperature  gradients  are  small.  Indeed,  for  a  mixing  layer  with  Mac  =  0.4,  the 
thermal  velocity  dilatation  is  much  smaller  than  the  acoustic  dilatation,  and  the  thermal  and 
acoustic  temperature  fluctuations  are  of  the  same  magnitude  (Figure  1).  Thus  in  this  case, 
the  thermal  component  of  the  flow  is  one  order  higher  than  that  for  the  boundary  layer.  This 
requires  minor  adjustment  to  the  decomposition  procedure. 


The  acoustic  and  nonacoustic  dilatation  and  temperature  fluctuations  were  determined  as 
in  the  boundary  layer  starting  from  the  equation 


P  dxj  Pr  Re  dxk  \  dx *  J 


(1) 


This  is  remarkably  successful,  even  at  quite  high  Mach  numbers.  For  example  at  Mac  —  0.8, 


Figure  2:  Contours  of  velocity  divergence  in  a  Ma  =  0.8  mixing  layer.  Shown  are  a)  the 
contribution  from  thermal  fluctuations,  and  b)  the  contribution  from  acoustic  fluctuations. 
The  shock  is  clearly  visible  in  b). 


where  in  the  two-dimensional  mixing  layer,  eddy  shocklets  appear,  this  decomposition  cor¬ 
rectly  isolates  the  shock  dilatation  in  the  acoustic  field  (Figure  2). 

Numerical  simulations  of  three-dimensional  compressible  shear  layers  are  now  under  way. 
These  simulations  are  begun  with  initial  conditions  taken  from  the  incompressible  simula¬ 
tions  of  Rogers  &  Moser  (1994).  In  addition  to  providing  the  turbulence  fields  required  to 
apply  the  decomposition  analysis.  These  new  simulations  allow  compressibility  effects  to 
be  evaluated  directly  by  comparison  to  the  evolutions  of  the  incompressible  fields. 


3  Other  Information 


Supported  under  this  grant  have  been  Stanislov  Borodai  and  Wai-Yip  Kwok,  both  graduate 
students  in  the  TAM  department  at  University  of  Illinois.  The  work  supported  here  was 
presented  at  the  1997,  1998  and  1999  meetings  of  the  Division  of  Fluid  Dynamics  of  the 
American  Physical  Society,  and  three  papers  have  resulted  from  this  work  (Borodai  & 
Moser,  2000,  Kwok,  Moser  &  Jimenez,  2000,  Guarini  et  al.,  2000). 
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(Received  20  March  1999  and  in  revised  form  21  January  2000) 

A  direct  numerical  simulation  of  a  supersonic  turbulent  boundary  layer  has  been 
performed.  We  take  advantage  of  a  technique  developed  by  Spalart  for  incompressible 
flow.  In  this  technique,  it  is  assumed  that  the  boundary  layer  grows  so  slowly 
in  the  streamwise  direction  that  the  turbulence  can  be  treated  as  approximately 
homogeneous  in  this  direction.  The  slow  growth  is  accounted  for  by  a  coordinate 
transformation  and  a  multiple-scale  analysis.  The  result  is  a  modified  system  of 
equations,  in  which  the  flow  is  homogeneous  in  both  the  streamwise  and  spanwise 
directions,  and  which  represents  the  state  of  the  boundary  layer  at  a  given  streamwise 
location.  The  equations  are  solved  using  a  mixed  Fourier  and  B-spline  Galerkin 
method. 

Results  are  presented  for  a  case  having  an  adiabatic  wall,  a  Mach  number  of 
M  =  2.5,  and  a  Reynolds  number,  based  on  momentum  integral  thickness  and  wall 
viscosity,  of  Ree>  =  849.  The  Reynolds  number  based  on  momentum  integral  thickness 
and  free-stream  viscosity  is  Ree  =  1577.  The  results  indicate  that  the  Van  Driest 
transformed  velocity  satisfies  the  incompressible  scalings  and  a  small  logarithmic 
region  is  obtained.  Both  turbulence  intensities  and  the  Reynolds  shear  stress  compare 
well  with  the  incompressible  simulations  of  Spalart  when  scaled  by  mean  density. 
Pressure  fluctuations  are  higher  than  in  incompressible  flow.  Morkovin’s  prediction 
that  streamwise  velocity  and  temperature  fluctuations  should  be  anti-correlated,  which 
happens  to  be  supported  by  compressible  experiments,  does  not  hold  in  the  simulation. 
Instead,  a  relationship  is  found  between  the  rates  of  turbulent  heat  and  momentum 
transfer.  The  turbulent  kinetic  energy  budget  is  computed  and  compared  with  the 
budgets  from  Spalart’s  incompressible  simulations. 


1.  Introduction 

The  study  of  supersonic  turbulent  boundary  layers  has  primarily  consisted  of 
experimental  investigations  with  a  few  recent  attempts  at  numerical  simulation.  The 
experimental  measurements  are  limited  to  basic  turbulence  quantities  and  by  the 
spatial  resolution  near  the  wall,  among  other  difficulties.  The  simulations  have  been 
hampered  by  large  cost  and  low  Reynolds  number.  The  goal  of  the  present  work 
was  to  identify  similarities  and  differences  between  compressible  and  incompressible 
boundary  layers,  as  well  as  to  test  the  applicability  of  Morkovin’s  hypothesis  and  the 
strong  Reynolds  analogy. 
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S.  E.  Guarini,  R.  D .  Moser,  K.  Shariff  and  A.  Wray 

l.L  Background  and  motivations 
Though  most  flows  encountered  in  nature  and  in  aerospace  applications  are  turbulent 
or  partially  so,  turbulence  remains  one  of  the  most  elusive  subjects  in  aeronautics. 
There  is  no  general  turbulence  theory  or  model.  With  the  addition  of  compress¬ 
ibility  the  turbulence  problem  becomes  even  more  complex.  For  instance,  in  bound¬ 
ary  layers  at  high  Mach  numbers,  large  temperature  gradients  develop  between  the 
wall  region  and  the  outer  layer.  These  gradients  result  in  large  variations  of  mean 
fluid  properties,  such  as  viscosity,  which  result  in  significant  changes  in  the  Reynolds 
number  across  the  wall-normal  direction.  As  the  Mach  number  becomes  hyper¬ 
sonic,  shocks  may  form,  dilatation  becomes  important,  and  baroclinic  terms  may  be 
significant. 

To  account  for  the  effects  of  compressibility,  many  theories  have  been  developed 
based  on  a  weakly  compressible  hypothesis  (see  the  review  by  Spina,  Smits  &  Robin¬ 
son  1994).  The  hypothesis  is  that,  at  moderate  free-stream  Mach  numbers  (M  <>  5 
according  to  Morkovin  1962),  dilatation  is  small  and  any  differences  from  incom¬ 
pressible  turbulence  can  be  accounted  for  by  fluid  property  variations  across  the 
layer.  This  has  long  been  appreciated  and  is  the  basis  of  the  Van  Driest  (1951,  1956) 
transformation  and  the  Morkovin  (1962)  hypothesis.  Morkovin  postulates  that  in 
the  weak  compressibility  regime  normal  stresses  will  obey  the  incompressible  scaling 
when  they  are  multiplied  by  the  local  mean  density  divided  by  the  free-stream  value. 
Even  at  moderate  free-stream  Mach  numbers  the  fluctuating  and  turbulence  Mach 
numbers  are  small  and  one  would  not  expect  eddy  shocklets  to  be  a  predominant 
feature  of  the  flow  field.  However,  at  higher  free-stream  Mach  numbers  the  turbulent 
velocity  fluctuations  are  more  likely  to  be  supersonic  leading  to  increased  compress¬ 
ibility  effects.  Further  support  for  the  weak  compressibility  assumption  is  provided  by 
the  fact  that  large-scale  structures  are  converted  at  0.9l/oo  (Spina,  Donovan  &  Smits 
1991)  which  results  in  a  small  relative  Mach  number  between  the  large  eddies  and  the 
mean  flow.  Finally,  another  measure  of  compressibility,  the  gradient  Mach  number 
(Sarkar  1995),  is  also  small  in  the  boundary  layer.  The  gradient  Mach  number  is 
based  on  the  velocity  difference  across  the  scale  of  an  eddy. 

The  validity  of  the  weak  compressibility  theories  in  compressible  boundary  layers 
has  been  checked  in  a  variety  of  experiments  over  the  years  (see  these  and  the 
references  therein:  Gaviglio  1987;  Smith  &  Smits  1993;  and  Elena  &  Gaviglio 
1993).  However,  in  most  experiments  the  data  reported  have  been  limited  to  simple 
turbulence  quantities  such  as  the  mean  and  RMS  velocity  and  temperature.  Detailed 
correlation  statistics  needed  to  directly  check  the  validity  of  Morkovin’s  hypothesis  are 
only  available  from  a  few  experiments.  The  data  from  the  direct  numerical  simulations 
reported  here  provide  an  opportunity  to  evaluate  these  theories  in  more  detail  than 
has  been  previously  possible. 

1.2.  Previous  simulations 

Although  there  have  been  numerous  experimental  investigations  of  the  compressible 
turbulent  boundary  layer,  there  have  been  relatively  few  attempts  at  direct  numerical 
simulation  of  this  flow.  To  date  there  have  been  three  such  attempts  known  to  us :  (i) 
Guo  &  Adams  (1994)  and  Adams  et  al  (1998);  (ii)  Rai,  Gatski  &  Erlebacher  (1995); 
and  (iii)  Hatay  &  Biringen  (1995). 

The  simulations  of  Guo  &  Adams  and  Adams  et  al  had  an  isothermal  wall  at  the 
laminar  adiabatic  wall  temperature  and,  like  the  current  simulations,  used  a  method 
that  transformed  the  spatially  evolving  boundary  layer  into  a  parallel,  streamwise 
homogeneous  flow.  To  obtain  their  transformed  parallel  shear  flow,  Guo  &  Adams 
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(1994)  require  that  the  spatial  mean  of  the  periodic  simulation  obey  the  parabolized 
time-mean  equations.  These  equations  contain  streamwise  derivatives  of  Reynolds 
stresses  which  are  obtained  by  performing  simulations  at  different  downstream  sta¬ 
tions.  This  approach  leads  to  forcing  terms  in  the  mean  equations  that  are  similar 
to  those  produced  by  Spalart’s  approach,  except  that  terms  involving  the  coordinate 
transformation  are  not  present. 

Guo  &  Adams  (1994)  and  Adams  et  al  (1998)  simulate  boundary  layers  at  three 
Mach  numbers  M  =  3, 4.5  and  6.  In  Guo  &  Adams,  these  simulations  were  performed 
in  relatively  small  spatial  domains;  in  particular,  for  M  =  3  the  domain  was  16  times 
smaller  in  streamwise  by  spanwise  area  than  those  reported  here.  The  effects  of 
these  small  domains  is  evident  in  their  statistical  results,  most  notably  the  two-point 
correlations.  Such  small  domains  introduce  considerable  uncertainty  regarding  the 
impact  of  the  domain  size  on  the  dynamics  of  the  boundary  layer,  though  experience 
with  the  minimal  channel  (Jimenez  &  Moin  1991)  suggests  that  basic  statistical 
quantities  such  as  turbulence  intensities  should  be  reasonably  accurate.  One  of  the 
major  design  considerations  of  the  current  simulations  was  to  avoid  the  uncertainties 
associated  with  such  small  simulation  domains. 

Unlike  the  current  simulations  and  those  of  Guo  &  Adams,  Rai  et  al  (1995) 
simulated  a  true  spatially  evolving  boundary  layer.  They  simulated  a  very  long 
streamwise  domain  from  a  laminar  inlet,  through  transition  to  a  fully  turbulent 
boundary  layer.  This  clearly  avoids  any  uncertainties  that  might  be  associated  with 
the  approximate  spatial  growth  treatment  used  here  and  in  Guo  &  Adams,  but 
the  cost  is  that  the  spatial  domain  that  must  be  simulated  is  spectacularly  large. 
As  a  result,  even  with  the  extremely  large  problem  size  they  were  able  to  run  (17 
million  nodes),  Rai  et  aV s  resolution  was  over  a  factor  of  3  coarser  in  the  streamwise 
direction  (measured  in  wall  units)  than  the  current  simulations.  Such  coarse  resolution 
introduces  uncertainties  that  are  different  from  those  associated  with  the  approximate 
spatial  growth  treatment.  Simulations  like  that  reported  here,  which  can  be  run  with 
much  better  spatial  resolution,  are  useful  for  assessing  the  impact  of  the  coarse 
resolution  that  must  be  used  in  a  simulation  like  that  of  Rai  et  al 

Finally,  Hatay  &  Biringen  (1995)  performed  a  parallel-flow  boundary  layer  cal¬ 
culation  at  M  =  2.5.  However,  the  data  they  present  suggest  that  the  turbulence  is 
not  being  sustained.  Indeed,  the  turbulence  intensities  appear  to  drop  significantly 
over  the  course  of  the  simulation.  It  is  possible  that  the  Reynolds  number  of  their 
simulation  is  too  low  to  sustain  turbulence.  The  authors  quote  a  Reynolds  number 
based  on  displacement  thickness  of  Res •  «  1000,  which  corresponds  to  a  momentum- 
thickness  Reynolds  number  of  Reo'  »  140.  Fernholz  &  Finley  (1980)  call  flows  in  the 
Reynolds  number  range  of  300  <  Re&  ^  6000  transitional,  based  on  an  analysis  of 
mean  velocity  profiles.  While  their  definition  of  transitional  flows  does  not  preclude 
having  sustained  turbulence  below  Rey  «  300,  the  Reynolds  number  in  Hatay  & 
Biringen’s  simulation  is  well  below  the  lower  limit  of  this  range. 

Clearly,  the  direct  numerical  simulation  of  a  compressible  turbulent  boundary 
layer  is  a  difficult  undertaking,  in  which  various  compromises  must  be  made  to  make 
the  simulation  practical.  In  the  research  reported  here,  we  have  pursued  the  most 
reliable  compressible  boundary  layer  simulation  that  we  were  able  to  do  with  current 
computation  capabilities,  choosing  good  spatial  resolution  and  adequate  domain 
size  over  true  spatial  evolution.  Such  a  simulation  allows  a  much  more  detailed 
analysis  of  compressibility  effects  in  the  boundary  layer  than  has  been  previously 
possible.  The  current  simulation  will  be  described  and  analysed  in  what  follows, 
which  includes  simulation  details  (§2),  turbulence  statistics  (§3),  Reynolds  analogies 
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(§4),  the  turbulent  kinetic  energy  budget  (§5)  and  conclusions  (§6).  The  details  of  the 
approximate  spatial  growth  treatment  are  given  in  the  Appendix. 

2.  Simulation  details 

The  details  of  the  current  simulation  are  provided  in  this  section.  The  parameters 
used  in  the  compressible  turbulent  boundary  layer  simulation  are  compared  with 
the  incompressible  boundary  layer  simulations  of  Spalart  (1988);  the  incompressible 
channel  simulation  of  Kim,  Moin  &  Moser  (1987);  and  the  compressible  boundary 
layer  simulations  of  Guo  &  Adams  (1994)  and  Rai  et  al  (1995). 

2.1.  Simulation  method 

One  difficulty  in  performing  compressible  turbulent  boundary  layer  simulations  is 
that  the  streamwise  direction  is  inhomogeneous.  This  precludes  the  use  of  periodic 
boundary  conditions,  and  as  a  result  FFTs,  in  this  direction.  Furthermore  it  neces¬ 
sitates  use  of  a  long  entrance  length  for  the  flow  to  adjust  from  artificial  in-flow 
conditions. 

Several  techniques  have  been  developed  to  address  one  or  both  of  these  issues 
for  incompressible  flow  (Spalart  &  Leonard  1987;  Spalart  1988;  Spalart  &  Watmuff 
1993;  Bertolotti,  Herbert  &  Spalart  (1992);  and  Lund,  Wu  &  Squires  1998).  Of  these 
we  utilize  the  one  developed  in  Spalart  &  Leonard  (1987)  and  Spalart  (1988).  Spalart 
recognized  that  the  slow  growth  of  the  boundary  layer  in  the  streamwise  direction 
makes  it  possible  to  treat  the  turbulence  as  approximately  homogeneous  in  this 
direction.  The  slow  growth  is  taken  into  account  by  using  a  coordinate  transformation 
and  a  transformation  of  dependent  variables  as  in  multi-scale  asymptotics.  The  result 
is  a  modified  system  of  equations  (Navier-Stokes  plus  some  extra  terms,  which  we  shall 
call  ‘slow  growth  terms’)  that  is  homogeneous  in  both  the  streamwise  and  spanwise 
directions,  and  which  represents  the  state  of  the  boundary  layer  at  a  given  streamwise 
location  (or,  equivalently,  a  given  thickness).  Using  Spalart’s  method,  the  boundary 
layer  can  be  simulated  separately  at  each  streamwise  station.  A  detailed  description 
of  this  method  and  the  modified  set  of  equations  can  be  found  in  the  Appendix. 

The  resulting  equations  are  solved  using  a  mixed  Fourier-spectral  and  B-spline- 
Galerkin  method  (Guarini  1998).  The  dependent  variables  (specific  volume,  a  = 
1/p;  momentum,  m  ~  pu;  and  pressure,  p )  are  expanded  in  terms  of  a  Fourier 
representation  in  the  horizontal  directions  and  a  third-order  (quadratic)  B-spline 
representation  in  the  wall-normal  direction.  The  Fourier  directions  are  de-aliased 
using  the  3/2-rule,  where  all  nonlinear  terms  are  calculated  using  3/2  times  the 
number  of  modes  used  to  advance  the  solution.  Quadratic  nonlinearities  are  fully  de- 
aliased  using  this  rule  while  higher-order  nonlinearities  are  only  partially  de-aliased. 
B-splines  have  a  variety  of  good  numerical  properties,  and  have  been  used  successfully 
in  the  incompressible  pipe  flow  simulation  of  Loulou  (1996)  and  the  compressible  jet 
of  Rao  (1997).  B-splines  have  high  resolving  power,  allow  easy  implementation  of 
boundary  conditions,  and  allow  the  use  of  stretched  grids.  More  details  on  B-splines 
may  be  found  in:  De  Boor  (1978),  Kravchenko,  Moin  &  Moser  (1996),  and  Shariff* 
&  Moser  (1998).  Their  use  in  the  present  work  is  described  in  Guarini  (1998).  In 
the  wall-normal  direction,  Giles’  (1989,  1990)  second-order  non-reflecting  boundary 
conditions  are  used  at  the  free-stream  boundary  and  adiabatic  no-slip  boundary 
conditions  are  used  at  the  wall.  This  combination  of  splines  and  Fourier  methods 
produces  a  very  accurate  numerical  method.  For  the  time  discretization  the  mixed 
implicit /explicit  method  of  Spalart,  Moser  &  Rogers  (1991)  is  used.  All  terms  are 
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Sim. 

M 

Ree 

Nx  X  N; 

Li  X  Lt 

Ax+  x  A z+ 

SI 

0.0 

225 

— 

_ 

— 

S2 

0.0 

300 

85x64 

2680+  x  670+ 

31.5+  x  10.5+ 

S3 

0.0 

670 

171  x  128 

4900+  x  1225+ 

28.7+  x  9.50+ 

S4 

0.0 

1410 

288  x  213 

11400+  x  2850+ 

39.6+  x  13.4+ 

KMM 

0.0 

— 

192  x  160 

2300+  x  1150+ 

12.0+  x  7.00+ 

GA 

3.0 

3015 

192  x  144 

527+  x  300+ 

2.74+  x  2.08+ 

GA 

4.5 

2618 

180  x  144 

260+  x  155+ 

1.44+  x  1.08+ 

GA 

6.0 

2652 

180  x  128 

229+  x  137+ 

1.27+  x  1.07+ 

Rai 

2.25 

6000 

971  x  321 

full  spatial 

27.0+  x  10.4+ 

Present 

2.5 

1577 

256  x  192 

2269+  x  1134+ 

8.86+  x  5.91+ 

Table  1.  Comparison  of  parameters  used  in  the  incompressible  simulations  of  Spalart  (1988)  (S1-S4) 
and  Kim  et  al  (1987)  (KMM);  the  compressible  simulations  of  Guo  &  Adams  (1994)  (GA)  and 
Rai  et  al.  (1995);  and  the  present  simulation. 


treated  explicitly  except  for  the  highest  wall-normal  derivative  in  the  viscous,  pressure 
gradient,  and  ‘acoustic  coupling  term’ 


yap 


dm2 

dx29 


(2.1) 


that  appears  in  the  pressure  equation.  In  the  implicit  treatment,  non-constant  coeffi¬ 
cients  that  vary  in  the  horizontal  directions  cannot  be  easily  treated.  Both  the  viscous 
and  acoustic  coupling  terms  are  split  into  a  term  with  coefficients  (viscosity  and  op 
respectively)  varying  in  the  wall-normal  direction  only,  which  is  treated  implicitly, 
and  the  remainder,  which  is  treated  explicitly. 


2.2.  Choice  of  parameters 

A  turbulent  boundary  layer  at  a  Mach  number  of  2.5  and  a  Reynolds  number  based 
on  displacement  thickness  of  Res •  =  6258  was  simulated.  This  results  in  a  Reynolds 
number  based  on  momentum  integral  thickness  and  wall  viscosity  of  Ree '  «  849.  The 
Reynolds  number  based  on  momentum  integral  thickness  and  free-stream  viscosity 
was  Ree  «  1577.  The  Mach  number  was  chosen  because  of  the  availability  of 
experimental  data  and  because  it  is  in  a  range  where  we  might  begin  to  see  some 
compressibility  effects. 

There  are  two  important  sets  of  parameters,  the  grid  size  ( Nx  x  Ny  x  Nz)  and 
the  domain  size  (Lx  x  Ly  x  Lz\  that  determine  the  overall  quality/accuracy  of  the 
simulation.  The  coordinate  system  is  oriented  such  that  the  x-,  y-,  and  z-directions 
are  the  streamwise,  wall-normal,  and  spanwise  directions,  respectively  (in  this  paper 
we  use  a  mixed  index  or  symbol  notation  where  xj,  x2,  and  X3  correspond  to  x, 
y,  and  z,  respectively).  The  current  simulation  has  256  x  209  x  192  ( Nx  x  Ny  x  Nz) 
Fourier-B-spline  modes  and  a  domain  size  of  2269+  x  875+  x  1134+  ( Lx  x  Ly  x  Lz), 
where  y+  =  yuT/vw.  Here  vw  is  the  kinematic  viscosity  at  the  wall  and  uT  is  the 
friction  velocity  (t w/Pw)1/2,  where  tw  and  pw  are  the  shear  stress  and  density  at  the 
wall,  respectively.  The  domain  and  grid  parameters  were  selected  to  provide  sufficient 
resolution  in  a  domain  that  is  large  enough  to  eliminate  most  finite  domain  size 
effects.  That  this  is  the  case  is  demonstrated  below. 

One  way  to  assess  the  adequacy  of  the  resolution  and  domain  size  is  by  comparison 
to  DNS  of  similar  flows.  The  resolution  and  domain  size  used  in  the  current  simulation 
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Figure  1.  One  dimensional  energy  spectra  £u#Ua.  Plotted  versus  (a)  kx  at  y+  =  4;  ( b )  kx  at  y+  =  80; 

(c)  kz  at  y+  =  4;  (d)  kz  at  y+  =  80. - ,  Streamwise  velocity  component; - ,  wall-normal 

velocity  component; - ,  spanwise  velocity  component. 


are  compared  with  the  incompressible  boundary  layer  simulations  of  Spalart  (1988) 
(S1-S4),  the  incompressible  channel  flow  of  Kim  et  al  (1987)  (KMM)  and  the 
compressible  boundary  layer  simulations  of  Guo  &  Adams  (1994)  (GA)  (table  1). 
The  resolution  required  in  the  M  =  3  simulation  of  Guo  &  Adams  was  significantly 
finer  than  that  of  the  incompressible  simulations  or  the  present  simulation.  The 
reason  for  Guo  &  Adams’  extremely  fine  resolution  is  not  clear.  However,  the  need 
for  increased  resolution  of  the  current  simulation  relative  to  KMM  is  due  to  sharp 
density  gradients  present  in  the  compressible  flow.  KMM  is  used  for  comparison 
because  their  resolution  is  better  than  Spalart’s  simulations,  as  determined  by  the 
drop  in  energy  spectra.  The  adequacy  of  the  spatial  resolution  was  confirmed  by 
examining  the  spectra,  where  E„aUa  is  the  energy  spectrum  for  the  velocity  component 
Ui.  Examples  are  shown  in  figure  1.  These  spectra  and  those  at  other  y-locations  suggest 
the  resolution  is  adequate.  Another  indication  of  the  adequacy  of  the  resolution  is 
the  value  of  kmaxri,  where  kmax  is  the  maximum  wavenumber  in  x  and  rj  is  the 
local  Kolmogorov  scale.  The  maximum  and  minimum  of  this  value  in  the  current 
simulation  are  1.6  and  0.5,  respectively,  which  is  considered  adequate.  For  comparison 
the  simulation  of  Kim  et  al  (1987)  had  values  of  1  and  0.4  for  the  maximum  and 
minimum,  respectively. 


Figure  2.  Two-point  correlations  QUaUa  -  Plotted  versus  (a)  Sx/S *  at  y+  =  4;  (6)  <5jc/<5*  at  y+  =  150;  (c) 

<52/<5*  at  y+  =  4;  (d)  <52/<5*  at  y+  =  150. - ,  Streamwise  velocity  component; - ,  wall-normal 

velocity  component; - ,  spanwise  velocity  component. 


The  periodic  domain  size  of  the  current  simulation  was  selected  to  ensure  that  the 
streamwise  and  spanwise  two-point  correlations  are  nearly  zero  for  large  separations, 
where  the  two-point  correlation  for  the  velocity  component  ux  is  QUaUa-  As  is  evident 
in  figure  2,  the  near-wall  correlations  are  indeed  near  zero  for  large  separations, 
though  they  could  be  better  for  the  streamwise  component.  However,  far  from  the 
wall  (y+  =  150)  low-level  large-scale  coherence  is  evident  in  the  correlation,  perhaps 
due  to  acoustics  as  suggested  by  Coleman,  Kim  &  Moser  (1995). 

In  the  wall-normal  direction,  the  B-splines  were  defined  on  a  non-uniform  set  of 
knots  (grid  points)  which  is  given  by 

t>  =  (Aymax  -  A ym„)  -  l)/o|  +  Aym^Nm,  (2.2) 

where  a  =  0.14,  Nm  —  1.09,  Nk  —  207,  A ymin  =  1.0,  and  A ymax  =  110.0.  With  this 
distribution  there  were  13  grid  points  in  the  first  9  wall  units,  including  the  grid  point 
at  the  wall.  The  minimum  grid  spacing  in  the  wall-normal  direction  was  0.48+  wall 
units  while  the  maximum  was  7.8+  wall  units  at  the  free-stream  boundary,  which  was 
located  at  y+  =  875. 
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Figure  3.  Turbulent  and  fluctuating  Mach  numbers  as  functions  of  y/S: - ,  turbulent  Mach 

number,  Mt; - ,  RMS  of  Mach  number  fluctuations,  (M'2)l/2; - ,  turbulent  Mach  number 

with  spanwise  component  removed. 


2.3.  A  note  on  averaging 

In  the  results  that  follow,  both  Reynolds  and  Favre  averaging  are  used  depending 
on  simplicity  of  presentation  and  conventions  used  in  the  papers  to  which  we  are 
comparing.  In  each  case  care  will  be  taken  to  distinguish  between  the  two. 

The  Reynolds  average  of  /  over  the  x-  and  z-directions  will  be  denoted  by  /,  and 
fluctuations  about  this  mean  will  be  denoted  by  The  Favre  average  over  the  x- 
and  z-directions,  /,  is  a  density-weighted  average: 

]=^r.  (2-3) 

P 

Fluctuations  about  the  Favre  average  will  be  denoted  by  /". 


2.4.  The  turbulent  Mach  number 

One  convenient  measure  of  the  turbulent  compressibility  effects  is  the  fluctuating 
Mach  number,  M',  which  is  just  the  RMS  fluctuation  of  the  Mach  number.  A  similar 
quantity  is  the  turbulent  Mach  number  given  by 


Mt  = 


a 


(2.4) 


Morkovin  (1962)  suggests  that  the  turbulence  is  only  weakly  affected  by  compressibil¬ 
ity  provided  M'  <;  0.2  (0.3  according  to  Spina  et  al  1994).  Despite  the  relatively  low 
Mach  number,  the  peak  values  of  M'  and  Mt  in  the  simulation  are  approximately  0.3 
(figure  3).  Nevertheless,  as  will  be  shown  in  §  3,  Morkovin’s  density  scaling  and  the  Van 
Driest  transformation  still  apply.  The  shoulder  in  figure  3  is  at  (M/2)1/2  «  0.25,  while 
the  experimental  data  shown  in  Spina  et  al  (1994),  for  =  2.3  and  Re0  =  4700, 
have  a  shoulder  at  (Af'2)1/2  «  0.20.  The  reason  for  the  difference  is  most  likely  the 
neglect,  in  the  experiments,  of  the  spanwise  velocity  in  calculating  the  Mach  num¬ 
ber.  As  one  can  see  in  figure  3,  the  turbulent  Mach  number  is  consistent  with  the 
experimental  value  of  0.2  when  the  spanwise  velocity  component  is  not  included. 

A  measure  of  intrinsic  compressibility  is  the  ratio  of  mean-square  dilatation  fluc¬ 
tuations  to  mean-square  vorticity  fluctuations: 


(du'i/dxMdu'j/dx^/w1^. 


(2.5) 
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Figure  4.  Comparison  of  the  computational  C/  with  the  experimental  data  of  Coles  (1954):  *, 

simulation  data  point;  o,  experimental  data  points  (2.2  ^  M  <  2.8); - ,  Van  Driest  II  (Bardina 

et  al  1997). 


This  ratio  measures  the  level  of  compressibility,  as  given  by  dilatation,  relative  to  the 
turbulent  motion,  as  given  by  enstrophy.  In  the  simulation  we  find  that  this  ratio  is 
approximately  5  x  10-4  throughout  the  boundary  layer. 


3.  Turbulence  statistics 

In  this  section,  several  turbulence  statistics  are  examined  to  evaluate  their  consis¬ 
tency  with  accepted  experimental  and  computational  results. 

To  obtain  statistics,  averages  are  computed  over  the  streamwise  and  spanwise 
directions  of  each  field;  then  an  ensemble  average  over  55  fields  spanning  31.33  time 
units  was  calculated.  Time  is  non-dimensionalized  by  S* /a^.  The  flow  was  determined 
to  be  stationary  when  several  quantities  (C/,  0,  wT,  Re'e,  and  Tw)  began  to  oscillate 
about  a  mean  value. 


3.1.  Mean  flow 

The  skin  friction  coefficient  is  defined  as 


In  the  simulation  the  skin  friction  coefficient  was  found  to  be  Cf  «  0.00282.  There 
are  very  few  experimental  studies  at  the  low  Reynolds  number  of  the  simulation. 
However,  the  simulation  compares  favourably  with  the  experimental  results  compiled 
by  Coles  (1954)  and  the  skin  friction  correlation  given  in  Bardina,  Huang  &  Coakley 
(1997)  based  on  the  Van  Driest  II  skin  friction  transformation  (figure  4). 

The  Van  Driest  transformed  velocity,  Uc,  is  plotted  in  wall  units  in  figure  5.  Uc  is 
defined  as 

Uc=  [  (TJT)l'2dU.  (3.2) 

Jo 

Experiments  have  shown  that  Uc  satisfies  the  same  scaling  laws  as  the  mean  streamwise 
velocity  in  incompressible  flow.  On  the  plot  we  have  included  the  linear  sub-layer 
relation,  17+  =  y+,  the  standard  log-law,  and  a  composite  profile  that  consists  of 
Reichardt’s  (1951)  inner  layer  profile  and  Finley’s  wake  function  (Cebeci  &  Bradshaw 
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Figure  5.  Van  Driest  transformed  velocity  in  wall  units: - ,  DNS  time  average; - ,  linear 

sub-layer  17+  =  y+; - ,  log-law  17+  =  (1/0.40)  In  (yuT/vw)  +  4.7; . ,  Reichardt’s  profile  with 

Finley’s  wake  function  ; - ,  Coles’  profile. 


1977): 


1  _  Q-yu t/foivw)  __  f  yUx  A  Q-yu tVK) 

\mvwJ 


+l  [(f)2  -  (f )’ + m  (f )!  -4n  (!)]  •  <33> 

where  <5  is  the  y  location  at  which  Uc(y  =  3)  =  UCn9  C\  —  — (l/jc)ln(fc)  +  C,  r\\  —  11, 
and  b  =  0.33.  This  profile  is  used  in  the  Appendix  to  calculate  the  slow-growth  terms. 
In  addition,  the  more  commonly  cited  Coles’  profile  (Coles  1956)  is  shown.  In  the 
region  of  30  <  y+  <  70,  the  simulation  data  fall  on  the  log-law  curve,  where  we 
have  chosen  the  constants  k  =  0.40  and  C  ~  4.7  for  the  plot;  k  was  determined 
by  finding  the  minimum  of  y+(dU+/dy+)  as  a  function  of  y+.  Using  this  value  and 
the  location  of  the  minimum,  C  was  then  calculated.  These  values  compare  well 
with  those  of  Spalart  (1988),  and  k  is  within  the  range  of  values  quoted  in  the 
literature  (see  Smits  &  Dussauge  1996).  At  low  Reynolds  numbers  the  log- region 
becomes  vanishingly  thin  making  the  determination  of  k  and  C  difficult.  There  is  also 
some  disagreement  as  to  whether  or  not  the  values  are  Reynolds  number  dependent 
at  the  low  Reynolds  number  of  the  simulation  (Spalart  1988).  The  value  of  77 
for  the  time-averaged  profile  was  77  =  0.25,  determined  from  the  equations  in  the 
Appendix.  Reichardt’s  ‘basic’  profile  with  Finley’s  wake  function  gives  a  rather  good 
representation  of  Uc  throughout  the  boundary  layer.  Since  this  is  the  profile  shape 
assumed  when  computing  the  slow  streamwise  derivatives  (see  the  Appendix),  the 
good  agreement  implies  that  the  assumed  profile  does  not  introduce  a  significant 
error  in  the  computation  of  these  derivatives.  It  is  interesting  to  note  that  if  we  use 
U  rather  than  the  transformed  velocity  Uc  the  values  of  k  and  C  are  0.477  and  2.64, 
respectively,  which  are  quite  different  from  the  incompressible  values. 


3.2.  RMS  velocity,  pressure,  and  vorticity 

When  normalized  by  ux  the  turbulence  intensities  from  the  current  compressible 
boundary  layer  are  lower  than  the  intensities  from  the  incompressible  boundary 
layer  (figure  6a).  Morkovin  (1962)  predicted  that  scaling  by  the  square  root  of 
the  mean  density  profile  should  collapse  RMS  data  for  the  streamwise  velocity 
component  and  possibly  the  spanwise  and  wall-normal  components.  When  this  scaling 
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Figure  6.  RMS  velocity  profiles  plotted  versus  y/5 ,  (a)  scaled  by  l/ur,  and  ( b )  scaled  by  y/(p/pw)/uT. 
Lines  correspond  to  the  compressible  simulation  and  symbols  are  used  for  Spalart’s  incompress¬ 
ible  simulations: - ,  streamwise  velocity  component; - ,  wall-normal  velocity  component; 

- ,  spanwise  velocity  component;  +,  Spalart  (Reg  =  1410);  *,  Spalart  (Ree  —  670). 


is  used  there  is  good  agreement  for  all  three  velocity  components  (figure  6b).  In 
the  experiments  good  collapse  is  obtained  for  the  streamwise  component  but  the 
experiments  are  inconclusive  with  respect  to  the  collapse  of  the  other  two  components 
(Smits  &  Dussauge  1996).  Smits  &  Dussauge  (1996)  attribute  this  to  both  the 
difficulty  in  measuring  the  spanwise  and  wall-normal  components  and  the  scarcity  of 
measurements  of  these  two  components. 

There  are  two  additional  points  which  need  to  be  mentioned  in  connection  with 
figure  6.  The  first  is  that  Ree ,  based  on  local  viscosity  in  the  compressible  simulation, 
varies  across  the  boundary  layer  from  849  to  1577.  Spalart’s  (1988)  two  simulations 
at  Ree  =  670  and  1410  span  this  range.  This  is  important  becaus#e  in  y/5  units, 
the  location  of  the  peaks  in  the  intensities  moves  toward  the  wall,  since  it  remains 
approximately  fixed  in  wall  units.  Further,  Spalart  showed  that  at  these  low  Reynolds 
numbers,  the  magnitude  of  the  peaks  in  intensities  increases  with  Reynolds  number. 
The  second  issue  concerns  the  choice  of  <5  in  figure  6.  As  discussed  by  Spalart,  the 
collapse  of  the  data  over  a  wide  range  of  Reynolds  numbers  is  sensitive  to  the  choice 
of  5.  In  making  our  comparison,  we  used  a  definition  of  5  based  on  the  composite 
profile  of  (3.3)  and  made  no  effort  to  find  a  definition  that  would  better  collapse  the 
data.  This  might  account  for  the  differences  that  are  evident  at  y  >  0.65<5. 

In  the  compressible  simulation  the  pressure  fluctuations,  when  scaled  by  pwu 2, 
are  larger  than  those  found  in  the  incompressible  simulations  through  most  of  the 
boundary  layer  (figure  7).  The  peak  pressure  fluctuations  are  larger  than  Spalart’s 
Ree  =  670  simulation  and  occur  at  nearly  the  same  location.  The  value  of  the  RMS 
pressure  at  the  wall  and  at  the  peak  are  (p'rms)w/(pwu2)  «  2.7  and  (p'rms)max/(pwu2)  «  3.0, 
respectively.  In  the  free  stream  the  pressure  fluctuations  approach  ( p'rms)oo/(pwu ?)  « 
0.47.  This  is  comparable  to  the  value  of  the  radiated  pressure  measured  by  Laufer 
(1964)ofO/JJM)*0.4. 

RMS  vorticity  profiles  for  the  present  simulation  agree  very  well  with  those  found 
in  Spalart’s  incompressible  simulations  when  normalized  by  u2/vw  (figure  8)  and 
plotted  in  wall  units,  with  the  compressible  results  being  slightly  larger  than  the 
incompressible  simulations  away  from  the  wall.  Note  that  from  Spalart’s  data  it  is 
clear  that  near  the  wall  the  maximum  may  be  Reynolds  number  dependent.  The 
near-wall  RMS  vorticity  is  shown  plotted  in  wall  units  here  because  the  wall  scaling 
is  known  to  (approximately)  collapse  such  profiles  for  different  Reynolds  numbers  in 
incompressible  flows,  and  indeed  in  y/5  units  our  data  did  not  collapse  with  that  of 
Spalart’s. 
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Figure  7.  RMS  pressure  fluctuations  versus  y/S : - ,  compressible  DNS;  +,  Spalart 

(Ree  =  1410);  *,  Spalart  (Ree  =  670);  x,  Spalart  (Ree  =  300). 


Figure  8.  RMS  of  vorticity  fluctuations.  Lines  correspond  to  the  compressible  simulation : 
- ,  G)f; - ,  g>22; - ,  GU32;  +,  Spalart  (Ree  =  1410);  *,  Spalart  (Reg  =  670). 


3.3.  Reynolds  shear  stress 
Reynolds  averaging  the  momentum  equation, 


dpui 

dt 


dp  1  dxjj_ 

dxi  Re  cxj  ’ 


we  get 


dpiii 


dp 


1  dzu 


-W—W,  + 


The  sum  of  the  Reynolds  shear  stress  and  mean  shear  stress  terms  is 


d__ 

dx2 


(3.4) 

(3.5) 


(3.6) 


In  incompressible  boundary  layers  a  constant  total  stress  region  is  observed  near  the 
wall.  The  constant-stress  region  is  consistent  with  the  law  of  the  wall  and  is  also 
present  in  the  current  simulation  (figure  9).  The  constant-stress  region  extends  to  30 
or  40  wall  units  above  the  wall. 

As  was  the  case  for  the  turbulence  intensities,  uV  agrees  with  the  incompressible 
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Figure  9.  Simulation  results  for  Reynolds  stress,  mean  shear  stress,  and  total  stress  versus  yux/vw: 

- ,  total  shear  stress; . ,  Reynolds  stress  (pw'/wj); - ,  mean  shear  stress  (] Ji(du/dy )); 

- ,  stress  correlations  (p'[(du'/dy)  +  ( dv'/dx )]). 


ylt 5  ylS 

Figure  10.  uV  versus  y/S ,  (a)  scaled  by  l/i£  and  (b)  scaled  by  (p/pw)/u*: - ,  compressible 

DNS;  +,  Spalart  ( Ree  =  1410);  *,  Spalart  ( Ree  =  670). 


results  when  scaled  by  the  local  mean  density  (figure  10).  This  agreement  is  not  as 
close  as  for  the  intensities,  but  this  is  primarily  due  to  the  square  root  in  the  definition 
of  the  intensities. 


4.  Reynolds  analogies 

For  incompressible  laminar  boundary  layers,  the  similarity  of  the  momentum  and 
energy  equations  allows  one  to  approximately  relate  quantities  pertaining  to  heat 
transfer  with  quantities  pertaining  to  momentum  transfer.  O.  Reynolds  discovered 
this  principle  in  its  simplest  form.  The  ‘Reynolds  analogy’  has  been  extended  with 
additional  approximations  to  the  compressible  and  turbulent  cases.  Morkovin  (1962) 
suggests  that  a  Reynolds  analogy  might  apply  to  compressible  turbulence,  a  concept 
known  as  the  ‘strong  Reynolds  analogy  (SRA)’.  More  recently,  other  expressions  of 
a  Reynolds  analogy  have  been  formulated  by  several  authors  (Gaviglio  1987  and 
Huang,  Coleman  &  Bradshaw  1995),  and  these  will  also  be  studied  below. 

4.1.  The  strong  Reynolds  analogy 

To  investigate  the  validity  of  the  strong  Reynolds  analogy,  and  its  consequences,  a 
brief  review  of  its  derivation  and  a  critical  examination  of  the  underlying  assumptions 
are  given  in  this  section.  The  analogy  is  based  on  the  observation  that  the  transport 
equations  for  mean  velocity  and  mean  total  enthalpy,  ht~CpT +uf/2,  for  a  stationary, 
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zero-pressure-gradient  boundary  layer, 

_ du\  _ du\  d  f  _ 


_dui 

dx2 


K'Mzj- 

dht  dht  d  (  p  dh,  — 

pU^  +  pU^2=^{rr372-pU^^ 


(4.1«) 

(4.1*0 


have  the  same  form  if  the  Prandtl  number  is  1.  If  the  Prandtl  numbers  are  different, 
the  equations  are  still  of  the  same  form  provided  the  molecular  diffusivity  can  be 
neglected,  which  is  true  in  a  turbulent  boundary  layer,  except  very  near  the  wall.  The 
two  equations,  however,  have  different  boundary  conditions. 

The  total  temperature,  T„  is  defined  by  the  relationship  h,  =  CpTt.  If  one  assumes 
that  gradients  of  mean  total  temperature  and  velocity  in  x  can  be  neglected  and 
further  (without  justification)  that 


T"  =  Cu", 


(4.2) 


one  can  eliminate  the  turbulent  terms  in  the  mean  transport  equations  and  solve  for 
ft  in  terms  of  iq, 

f,  =  Ctq  +  D.  (4.3) 

This  result  was  first  stated  in  a  slightly  different  form  by  Young  (1953)f  and  later 
by  Morkovin  (1962),  who  called  assumption  (4.2)  the  strong  Reynolds  analogy.  Since 
dft/dy  =  0  at  an  adiabatic  wall  and  the  velocity  gradient,  dui/dy,  is  non-zero  at  the 
wall,  it  follows  that  the  constant,  C,  must  be  zero.  This  implies  that  the  mean  total 
temperature  is  constant  with  value  flw  and  the  total  temperature  fluctuations  are 
zero.  In  the  current  simulation  the  maximum  deviation  of  the  mean  total  temperature 
from  a  constant  is  about  7%,  thus  approximately  verifying  the  result  for  the  mean. 
However,  in  the  simulation,  RMS  total  temperature  fluctuations  are  comparable  in 
magnitude  to  the  static  temperature  fluctuations  (see  figure  11),  and  are  thus  not 
negligible.  The  discussion  that  follows  addresses  the  validity  of  (4.2),  with  C  =  0,  and 
the  relationships  derived  from  it. 

The  fact  that  measured  total  temperature  fluctuations  are  not  negligible  was  rec¬ 
ognized  by  Morkovin  (1962).  Nonetheless,  relations  derived  assuming  that  they  are 
negligible  have  been  widely  used  in  the  literature.  These  relations  can  be  obtained  by 
writing  the  definition  of  total  temperature  and  subtracting  its  Favre  mean: 


Cp  T"  =  CVT"  +  u,<  + 


u"«|'  u"u" 
~2  2~' 


(4.4) 


By  retaining  only  the  terms  that  are  linear  in  the  fluctuations,  and  assuming  that 
tqu,  >  u2«2  and  tqu"  »  u3u3,  we  get 


CpT"  =  CpT"  +  uxxl[. 


(4.5) 


So  far,  the  approximations  made  are  excellent.  One  can  verify  this  by  considering 
the  correlation  coefficient  Ruj(r,"-r")  between  streamwise  velocity  fluctuations  and  the 
difference  between  total  and  static  temperature  fluctuations.  The  correlation  coefficient 
differs  from  unity  by  less  than  0.9%  for  y/d  >  0.05,  showing  that  (4.5)  is  very  accurate. 

In  Morkovin’s  analysis,  the  SRA  is  invoked  to  argue  that  the  total  temperature 


t  Young  calls  (4.3)  and  (4.2)  solutions  to  the  equations.  This  terminology  is  imprecise  and  attaches 
too  much  legitimacy  to  (4.3)  and  (4.2).  We  prefer  not  to  identify  the  two  as  solutions,  but  rather  to 
say  that  the  assumption  (4.2)  implies  (4.3). 
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Figure  11.  Comparison  of  RMS  total  and  static  temperature  fluctuations: - ,  (T!2)l/2/Tt; 

- ,  (T^/T; - ,  (T*)l'2/T. 


fluctuations  are  negligible  compared  to  the  static  temperature  fluctuations,  and  thus 
(4.5)  becomes 

CPT"  +  mh4«0,  (4.6) 

which  is  not  valid,  as  discussed  above  (figure  11).  Substituting  yR/(y  —  1)  for  Cp  and 
defining  the  Mach  number,  Ma : 


M]=  “L, 

*  yRT 

(4.7) 

(4.6)  can  be  rewritten: 

T"  i /' 

4-  =  -(y_l)M2^. 

T  u  i 

(4.8) 

The  (questionable)  relations,  (4.6)  and  (4.8),  between  T"  and  u'{ 
consequences  which  are  given  by 

have  several  statistical 

(fJay/2/f  ^  ^ 

(y  -  ~  ’ 

(4.9a) 

&  1, 

(4.9  b) 

pr  __  pu'lu'Wf/dy )  ..  , 
pxi’1Tn{du\ldy) 

(4.9  c) 

These  are  three  of  the  five  SR  A  relationships  Morkovin  (1962)  presented.  Equations 
(4.9a)  and  (4.92?)  are  direct  consequences  of  (4.6),  whereas  (4.9c)  is  obtained  by 
multiplying  (4.6)  by  pu averaging,  and  assuming  that 


which  comes  from  the  mean  total  temperature  equation  when  it  is  assumed  that  ft 
is  constant  and  that  nonlinear  fluctuating  terms  and  velocity  components  other  than 
mi  are  small. 

From  the  simulation  results  we  see  that  equation  (4.9a)  is  satisfied  for  y/S  <  0.6 
(figure  12).  However,  in  the  same  region,  (4.92?)  and  (4.9c)  are  not  satisfied  (see  figures 
13  and  14).  The  correlation  R^t”  is  approximately  0.6  through  most  of  the  boundary 
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y/S 

Figure  12.  Test  of  the  strong  Reynolds  analogy,  as  expressed  by  (4.9 a). 


Figure  13.  Correlation  coefficient,  -/fy r/,  versus  y/S.  (a)  Comparison  with  compressible  experi¬ 
ments:  - ,  time  average  of  DNS;  x,  Debieve  (M  —  2.3,  Reg  =  5650)  (Gaviglio  1987);  *,  Debieve 

(M  =  2.3,  Reg  =  5650)  (Gaviglio  1987);  +,  Dussauge  (M  =  1.7,  Ree  =  5700)  (Gaviglio  1987);  o, 
Smith  &  Smits  (1993)  (M  =  2.9,  Reg  =  77  600).  (b)  Comparison  with  incompressible  experiments: 

- ,  time  average  of  DNS;  x,  Fulachier  (Reg  =  5000)  (Gaviglio  1987);  *,  Rey  (Gaviglio  1987); 

+,  Subramanian  &  Antonia  (1981)  Reg  =  990;  o,  Subramanian  &  Antonia  (1981)  Reg  =  7100; 
- ,  incompressible  simulation  of  Bell  &  Ferziger  (1993). 


yIS 

Figure  14.  Turbulent  Prandtl  number,  Prt ,  versus  y/S: - ,  time  average  DNS  data;  x, 

incompressible  heated  wall  simulations  of  Bell  &  Ferziger  (1993). 

layer,  and  the  turbulent  Prandtl  number  is  about  0.7  except  near  the  wall  ( y/d  <  0.2) 
where  it  is  near  1.  The  current  results  for  the  velocity-temperature  correlation  are  in 
disagreement  with  the  available  data  from  compressible  boundary  layer  experiments 
(figure  13),  which  show  the  correlation  to  be  close  to  1  (approximately  0.9),  in 
agreement  with  predictions  of  the  SRA. 
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The  reason  for  this  disagreement  between  simulations  and  experiments  is  not  clear, 
but  is  of  great  interest.  The  speculations  of  Gaviglio  (1987)  that  the  large  correlation 
values  observed  in  experiments  are  due  to  acoustic  waves  suggests  the  possibility 
that  the  acoustics  are  somehow  different  (weaker)  in  the  simulations.  However,  tests 
for  numerical  artifacts  (e.g.  damping)  in  the  simulated  acoustics  did  not  reveal  any 
problems,  and  our  radiated  sound  pressure  magnitude  compares  well  with  that  of 
Laufer  (1964).  Another  interesting  observation  is  that  in  the  current  simulations 
both  the  velocity-temperature  correlations  and  the  turbulent  Prandtl  number  agree 
reasonably  well  with  data  from  incompressible  boundary  layers,  both  experimental 
and  computational  (see  figures  13 b  and  14).  This  might  be  expected,  given  the  weak 
compressibility  of  the  boundary  layer  at  this  Mach  number.  This  agreement  at  least 
makes  plausible  the  proposition  that  the  current  simulation  results  are  physical, 
despite  disagreement  with  compressible  experiments.  The  basis  for  comparison  with 
the  incompressible  heated  wall  in  figure  13  (b)  is  that  the  heated  wall  produces  a 
temperature  gradient  of  the  same  sign  as  that  obtained  in  the  current  simulation. 

Further  analysis  suggests  that  the  key  to  the  discrepancy  is  the  magnitude  of 
the  total  temperature  fluctuations.  In  the  experiments,  it  is  total  temperature  and 
momentum  that  are  measured  directly,  with  other  quantities  computed  using  several 
approximate  relations.  By  using  these  relations  in  the  simulation  data  it  was  found 
that  (a)  they  yield  accurate  values  for  the  derived  quantities;  and  ( b )  when  the  total 
temperature  fluctuations  in  the  simulations  are  reduced  by  a  factor  of  2,  the  relations 
used  in  the  experiments  yield  correlations  R^t"  consistent  with  the  experiments. 
The  total  temperature  fluctuations  in  the  experiments  are  indeed  approximately  a 
factor  of  2  lower  than  the  simulations  (Kistler  1959).  A  factor  of  2  error  in  the 
total  temperature  is  much  larger  than  the  uncertainties  in  either  the  experiments  or 
the  simulation,  so  there  is  clearly  something  wrong  with  one  or  the  other  (or  both). 
A  potentially  useful  approach  to  settling  this  question  would  be  to  use  a  physical 
model  of  the  experimental  probes  in  the  simulations,  and  process  the  resulting  signals 
as  is  done  in  the  experiments,  to  see  if  results  consistent  with  the  experiments  are 
recovered.  Until  this  issue  is  resolved,  it  would  be  wise  to  view  both  the  experimental 
and  computational  results  on  these  quantities  with  some  suspicion. 

Of  the  three  relations  shown,  only  the  RMS  relation  (4.9a)  is  very  nearly  satisfied. 
The  question  arises  as  to  how  this  relationship  can  be  satisfied  even  though  total 
temperature  fluctuations  are  of  the  same  order  as  temperature  fluctuations.  This 
success  of  (4.9a)  can  be  explained  by  rearranging  the  definition  of  total  temperature 
fluctuations  (4.5)  as  follows: 


T”2  +  V*  _  2T/T"  1x2,,4«i'2 

- hfi - —  =  (7  -  DX4  Jp 


(4.11) 


The  condition  that  must  be  satisfied  for  (4.9a)  to  be  valid  when  total  temperature 
fluctuations  are  not  neglected  is 


T,"2  -  2  T"  T" 

fi  ’ 


(4.12) 


which  the  simulation  data  confirms  (figure  15).  While  all  the  terms  on  the  left-hand 
side  of  (4.11)  are  of  the  same  order,  in  the  inner  portion  of  the  boundary  layer  the 
term  on  the  left-hand  side  of  (4.12)  is  nearly  a  factor  of  6  greater  than  the  term  on 
the  right.  So  the  success  of  (4.9a)  is  due  to  a  relationship  between  total  and  static 
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Figure  15.  Comparison  of  terms  in  (4.12): - ,  T,,2/f 2; - ,  T/'2/f2; . ,  — 2T/,T"/t’2; 

- ,  (T72~277r7)/f2. 


Figure  16.  Comparison  of  R^'t"  from  (4.15)  to  simulation  data: - ,  time  average  of  DNS  data; 

- ,  profile  predicted  by  (4.15); - ,  profile  predicted  by  the  modified  Reynolds  analogy  of 

Huang  et  al.  (1995). 


temperature  fluctuations  rather  than  the  assumption  of  negligible  total  temperature 
fluctuations. 

Gaviglio  (1987)  has  shown  that  fluctuations  in  total  and  static  temperature  can  be 
directly  related  to  the  correlation  coefficient  Ru”t"  if  (4.9a)  is  assumed  to  be  valid. 
The  RMS  of  the  linearized  definition  of  total  temperature,  (4.5),  is 

(T[12)112  =  (t"2  +  +  2^t(T"2)1/2(u"2),/2Ru»t"J  ,  (4.13) 

where  the  last  term  was  written  in  terms  of  R^t”-  Now  if  we  say  that  (4.9a), 


(T«2)l/2  =  £L(u"2)1/2, 


is  empirically  valid  then  (4.13)  becomes 


rrtf2 

Ru'!T "  =  '  ==  ““  1* 
1  2  T”1 


(4.14) 


(4-15) 


As  expected  in  the  case  of  negligible  total  temperature  fluctuations  this  equation 
reduces  to  Ru»t«  =  -1.  The  value  of  Ru"T"  predicted  by  (4.15)  agrees  well  with  the 
actual  values  (figure  16). 
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Figure  17.  Comparison  of  —R^r,  —R^r,  and  —Rp>r  versus  y/5: - ,  —R,^ r ; - >  —RP'r> 

- ,  — symbols  showing  the  experimental  results  for  —Ru>:  v  are  the  same  as  in  figure  13a. 


Since  for  incompressible  flow  —R^v  is  significantly  less  than  unity,  the  correlation 
coefficient  between  streamwise  momentum  fluctuations  and  temperature  fluctuations, 
— Rmfj /,  will  also  be  much  less  than  unity  for  incompressible  flow  (the  two  are  in  fact 
equal).  For  the  compressible  simulation,  however,  —Rm'j1  is  much  closer  to  1  (figure 
17).  This  is  because  —Rm\r  is  a  weighted  average  of  RU'[T>  and  RP’T>>  and  Rp>r  is  very 
close  to  1  (dashed  curve  in  figure  17)  since  pressure  fluctuations  may  be  assumed 
to  be  small  compared  to  density  and  temperature  fluctuations  in  the  equation  of 
state  (Lele  1994).  The  contribution  of  Rpt  to  the  weighted  average  pushes  —Rm^r 
towards  unity.  The  close  correlation  between  streamwise  momentum  and  temperature 
fluctuations  may  indicate  a  greater  similarity  between  the  transport  equations  for 
turbulent  momentum  and  heat  transport  than  in  the  incompressible  case,  which  may 
be  due  to  a  reduction  in  the  importance  of  the  pressure  gradient  term. 


4.2.  ' Modified  Reynolds  analogies' 

Both  Gaviglio  (1987)  and  Huang  et  al  (1995)  point  out  that  for  flows  with  non- 
adiabatic  boundaries  the  agreement  between  (4.9a)  and  measurements  is  poor.  Both 
authors  propose  new  relationships  between  temperature  and  velocity  fluctuations 
which  have  been  called  ‘modified’  Reynolds  analogies.  Since  the  agreement  of  this 
relation  with  the  current  adiabatic  wall  simulations  is  not  perfect  either,  we  now  assess 
these  modified  Reynolds  analogies.  The  modified  Reynolds  analogies  of  Rubesin 
(1990),  Gaviglio  (1987)  (GSRA),  and  Huang  et  al  (1995)  (HSRA)  all  have  the  form 


(7V/2)l/2/f _ 1 _ 

(y  -  1  )MH^)l/2/ui  *  c(*  ^  a(dft/df)Y 


(4.16) 


If  a  =  0  and  c  =  1  then  the  standard  SRA  is  obtained.  For  all  the  modified 
expressions  a  =  1.  Rubesin  used  c  =  1.34.  Gaviglio  and  Huang  et  al  use  the  mixing 
length  hypothesis  to  obtain  c  =  1.0  and  c  =  Prh  respectively.  The  difference  between 
the  derivation  of  Gaviglio  and  Huang  et  al  is  that  Gaviglio  assumes  that  the  mixing 
lengths  for  temperature  and  velocity  fluctuations  are  equal.  In  figure  18  the  ratio  of 
the  left-hand  side  of  (4.16)  to  the  right-hand  side  is  plotted  for  the  SRA,  GSRA,  and 
HSRA,  The  version  presented  by  Huang  et  al  (1995)  is  in  excellent  agreement  with 
the  data  throughout  the  boundary  layer.  The  value  of  RU^T"  that  is  predicted  by  the 
modified  analogy  of  Huang  et  al  can  be  derived  by  substituting  equation  (4.16),  with 
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Figure  18.  Plot  of  the  Strong  Reynolds  analogy,  and  the  modified  Reynolds  analogies  of  Gaviglio 
(1987)  and  Huang  et  al.  (1995): - ,  SRA; - ,  GSRA; - ,  HSRA; . ,  Rubesin. 


c  —  Prh  into  (4.13)  to  obtain 


[Tta/T>*- 1]  _  Pr,  (  STA 

[2Prt(l-dft/df)]  2  V  df)  <V 


(4.17) 


The  profile  for  RU»T»  predicted  by  (4.17)  agrees  well  with  the  simulation  data  (figure 
16). 

Why  does  the  expression  of  Huang  et  al  work  so  well?  Substituting  the  definition 
of  the  turbulent  Prandtl  number,  (4.9c),  and  the  derivative  of  the  total  temperature, 


into  (4.16),  with  c  • 


8ft 

8y 

■  Pr„  we  obtain 


df_ 

dy 


puf2  T"  _ 
(T""2)l/2 


Ml  du\ 

Cp  dy  ’ 

(4.18) 

pu'^u’l 

(4.19) 

(uf)1/2' 

This  relationship  expresses  an  analogy  between  the  rates  of  turbulent  heat  and 
momentum  transfer  normalized  by  the  property  that  is  transported.  We  may  divide 
through  by  the  RMS  of  the  wall-normal  velocity  fluctuations  to  obtain  a  relationship 
between  the  correlation  coefficients  R^t "  and  Ruy>: 


Ru$T"  =  —Ru”u'±, 


(4.20) 


where  it  is  assumed  that  the  correlations  pV2'u"  and  p'u'^T”  are  small  compared  to 
pu^u'l  and  pu'^T".  The  simulation  results  indicate  that  the  two  correlation  coefficients 
are  very  nearly  equal  throughout  the  boundary  layer  (figure  19). 


4.3.  Turbulent  Prandtl  number 

Morkovin  was  aware  that  total  temperature  fluctuations  are  not  negligible  compared 
to  temperature  fluctuations  and  stated  that  another  set  of  expressions  could  be 
developed  by  assuming  that  pu'^T"  is  much  smaller  than  pu'^T".  In  the  lower  half  of 
the  boundary  layer  (y/S  <  0.5)  this  is  a  good  assumption  (figure  20). 

Using  this  assumption,  an  expression  can  be  developed  for  the  turbulent  Prandtl 
number.  Multiplying  (4.5)  by  pw2  and  averaging  gives 

pupf  =  pu'2'T"  +  ^-pu'X- 

Op 


(4.21) 
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Figure  19.  Comparison  of  RqTtt  and  —Ru»u»  versus  y/S: - ,  Ru^Tffl - ,  — 


Figure  20.  Plot  of  the  ratio  pu'^T"  /  pu'^T"  versus  y/S. 


Neglecting  pu'2'rf"  relative  to  pu2T"  yields 

pu'ju'l  __  __ Cp 

pufT"  “  wi ’ 

Substituting  (4.18)  into  (4.22)  gives 

Pr  =  (gjVgy)  =  ££  ( £1  _ 

'  pu^T"  (dui/dy)  dy  \dy  dy  ) 


(4.22) 


(4.23) 


This  prediction  agrees  with  the  simulation  data  in  the  inner  portion  of  the  boundary 
layer,  but  as  the  boundary  layer  edge  is  approached  the  agreement  becomes  poor 
(figure  21). 


5.  Turbulent  kinetic  energy  budget 

For  the  benefit  of  those  formulating  turbulence  models,  the  budgets  for  both 
the  Reynolds  stresses  and  the  turbulent  kinetic  energy  have  been  calculated.  In 
this  section,  the  turbulent  kinetic  energy  budget  is  presented  and  compared  with 
the  incompressible  simulations  of  Spalart  (1988).  The  Reynolds  stress  budgets  are 
presented  in  Guarini  (1998).  Favre  averages  are  used  in  the  analysis  to  simplify  the 
resulting  equations. 
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Figure  21.  Turbulent  Prandtl  number  Prt  vs.  y/<5: - ,  time  average  DNS  data; 

- ,  Prt  from  (4.23). 


The  turbulent  kinetic  energy  is  defined  as: 


-k  =  i  K< 


(5.1) 


2  p  ’ 

and  the  turbulent  kinetic  energy  equation  is,  after  assuming  homogeneity  in  the  x- 
and  z-directions,  given  by 


j-(pic)  +  u2-^—(pk)  —  P  +  T  +  n  +  D  —  <j)+Vc. 
at  ox  2 


The  symbols  are  defined  as 


. 1 dlli 

1  d 


T— 


n,n,+ni,.A.m+/A 


r,  S  U"  , 

°  ~  Sx2  ReX‘2' 


(5.2) 

(5.3a) 
(5.3  b) 
(5.3c) 
(5.3  d) 


(j)  = 


Re  dxj  ’ 


(5-3e) 


Vc  =  -u‘ 


dx2  1  dxi 


dx2 


(5.3/) 


The  terms  in  (5.2)  can  be  interpreted  as  follows:  the  left-hand  side  is  the  substantial 
derivative  of  the  turbulent  kinetic  energy  along  a  mean  streamline;  P  is  the  rate  of 
generation  of  turbulent  kinetic  energy  by  mean  velocity  gradients;  T  is  turbulent 
transport;  II  are  the  pressure  terms  (pressure  diffusion  and  pressure  dilatation, 
respectively);  D  is  viscous  diffusion;  —cj)  is  viscous  dissipation  per  unit  volume;  and 
finally,  Vc  includes  the  terms  that  arise  when  the  density  is  not  constant.  The  first 
two  terms  of  Vc  are  due  to  the  difference  between  the  Favre  and  Reynolds  average 
and  the  last  term  is  the  production  term  due  to  dilatation.  The  pressure  dilatation 
as  well  as  the  dilatational  dissipation,  which  are  not  included  in  Vc,  are  also  due 
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Figure  22.  Turbulent  kinetic  energy  budget:  - ,  simulation  data; - ,  Spalart’s  data 

Ree  =  1410; - ,  Spalart’s  data  Ree  =  670.  Symbols  are:  C,  Convective  term;  P,  Generation 

term;  T,  Turbulent  transport;  17,  Pressure  terms;  D,  Viscous  transport;  —0,  Viscous  dissipation. 


to  non-constant  density.  In  the  literature  the  dissipation  per  unit  mass  is  commonly 
referred  to  as  e.  Here  we  use  the  dissipation  per  unit  volume,  which  we  denote  as 
0  for  clarity.  The  dissipation  and  pressure  terms  are  in  a  form  similar  to  that  given 
in  Huang  et  al  (1995).  The  compressible  turbulent  kinetic  energy  budget  agrees  with 
the  incompressible  results  of  Spalart  at  two  different  Reynolds  numbers  (figure  22). 
The  generation  term  is  larger  than  in  the  incompressible  simulation,  and  the  viscous 
transport  and  turbulent  transport  terms  are  also  greater  in  magnitude  than  in  the 
incompressible  simulation.  Vc  is  small  and  has  not  been  included  on  the  plot  for 
clarity.  Its  maximum  value  is  a  factor  of  25  smaller  than  that  of  the  generation  term. 

The  effects  of  compressibility  on  the  dissipation  have  been  of  interest  in  the 
literature,  especially  in  the  context  of  compressible  turbulence  models  (Zeman  1990 
and  Sarkar  et  al  1991).  To  study  dissipation  in  the  current  simulation,  consider  0 
which  can  be  expanded  as 


=  JLM  (Ml  ,  M _ h.  Ml)  +  (Ml  +  M _  1 S- 

Redxi  \dxi  dx\  3  ll  dxk  J  Re  dxi  \dxi  dxi  3  ll dxk  J 


pf  du'j  (8Ui_  d_U}_  _  2  duk\ 
Redxi\dxi  dx^  3  ll dxk  J  ’ 


(5.4) 


where  the  three  terms  in  this  expression  will  be  referred  to  as  <j>i,  </>2,  and  <fo, 
respectively.  The  fluctuations  about  the  Favre  average  in  (5.3e)  have  been  replaced 
using  the  identity 


Re  dxi  Re  dx\ ' 


The  terms,  </>2  and  $3,  which  involve  viscosity  fluctuations,  are  negligible  compared 
to  in  this  simulation. 

The  first  term,  </>i,  can  be  decomposed  into  parts  that  are  more  amenable  to 
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Figure  23.  Comparison  of  dissipation  terms: - ,  — 0i  =  —  (<j)s  +  <Pi  +  <t>d)\ - ,  solenoidal 

dissipation  (— <j)s); - ,  dissipation  due  to  inhomogeneity  (—<£,);  . ,  dissipation  due  to 

dilatation  +,  Spalart’s  data  Ree  =  1410;  *,  Spalart’s  data  Reo  =  670. 


comparison  with  incompressible  flows  by  expressing  velocity  gradients  in  terms  of  the 
rate-of-deformation  tensor,  S'u,  and  spin  tensor,  Q'u,  which  is  related  to  the  vorticity 
through  Q'ijQ'tj  =  co'co[/2.  Simplifying  gives  for  (j)\ : 


0 


c ofa'i  +  2 


Re 


(  &  —n  o  \  4  du* 

I  "Z r W.'Mi —  2- Wl  — ]  -|-  —  — —  — , 

l  dxidxi  1  dxi  dxij  3  Redxtdxk 


(5.6) 


where  the  first  term  on  the  right  is  the  homogeneous  incompressible  dissipation,  or 
the  solenoidal  part  of  the  dissipation,  (j)s ,  the  second  term,  <j>h  is  due  to  inhomogeneity, 
and  the  third  term,  cj)d,  is  due  to  dilatation.  Both  the  dissipation  due  to  dilatation 
and  inhomogeneity  are  very  small  compared  to  the  solenoidal  dissipation  (figure  23). 
At  the  wall  the  dissipation  due  to  inhomogeneity  provides  a  very  slight  contribution 
to  the  total  dissipation.  The  compressible  result  agrees  well  with  the  incompressible 
results  of  Spalart. 

Finally,  we  consider  the  pressure  terms.  There  are  three:  pressure  diffusion  (77f), 
pressure  dilatation  (77^),  and  compressibility  (J7C).  The  pressure  dilatation  term  is  also 
associated  with  compressibility  effects,  since  the  dilatation  is  zero  for  an  incompressible 
flow.  The  three  pressure  terms  are 


nt  =  - 


di /2'p' 

HxT9 


nd  =  pf 


<M[ 

dxC 


nc  =  -hA 

OX  2 


(5.1  a-c) 


Of  the  three,  nd  and  nc  are  much  smaller  than  FIt  near  the  wall  (figure  24).  In 
fact,  the  sum  of  the  pressure  terms  is  almost  indistinguishable  from  FJt  in  the  wall 
region,  which  shows  that  the  compressibility  terms  have  very  little  effect  on  the  overall 
contribution  of  the  pressure  terms  to  the  turbulent  kinetic  energy  budget.  Away  from 
the  wall,  all  the  terms  are  small  and  the  compressibility  term  contributes  to  the  sum 
of  the  terms.  The  pressure  diffusion  term  is  larger  than  the  value  obtained  by  Spalart 
in  both  his  Re0  =  1410  and  Ree  ~  670  cases. 

The  results  and  analysis  given  here  show  that  at  M  =  2.5,  the  effects  of  compress¬ 
ibility  on  the  turbulent  kinetic  energy  balance  are  not  due  to  the  new  compressibility 
terms  that  appear  in  the  equations.  Instead,  the  effects  are  more  subtle,  quantitatively 
affecting  the  terms  that  appear  in  the  incompressible  case. 
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Figure  24.  Comparison  of  pressure  terms: - ,  +  /7d  +  I7C; - ,  77* ; - ,  nd\ 

. ,  77c;  +,  Spalart’s  data  Re$  =  1410;  *,  Spalart’s  data  Reg  —  670. 


6.  Conclusions 

A  direct  numerical  simulation  of  a  Mach  2.5  turbulent  boundary  layer  was  carried 
out  using  the  method  described  in  the  Appendix.  The  Reynolds  number  of  the 
simulation  was  Ree>  =  849.  Comparison  with  available  experiments  and  with  other 
simulations  (i.e.  those  of  Spalart)  suggest  that  the  current  simulation  provides  an 
accurate  description  of  a  compressible  turbulent  boundary  layer. 

It  was  shown  that  many  of  the  scaling  relations  used  to  express  compressible 
boundary  layer  statistics  in  terms  of  those  for  incompressible  boundary  layers  are 
consistent  with  the  current  simulation.  In  particular,  we  have  shown  that  the  Van 
Driest  transformed  velocity  behaves  much  the  same  as  the  streamwise  velocity  in 
the  incompressible  case.  There  was  a  small  logarithmic  region  with  k  ~  0.40  and 
C  =  4.7.  It  was  also  shown  that  the  RMS  velocity  fluctuations  are  collapsed  with 
incompressible  results  by  the  mean  density  scaling  suggested  by  Morkovin.  When  this 
scaling  is  applied,  the  data  from  the  current  simulation  agree  remarkably  well  with 
Spalart’s  Ree  =  670  and  Ree  =  1410  simulations.  The  mean  density  scaling  of  u'v' 
also  results  in  a  fairly  good  collapse  with  incompressible  results. 

An  inconsistency  with  the  standard  analysis  of  compressible  turbulent  boundary 
layers  was  found  in  that  the  total  temperature  fluctuations  were  of  the  same  order 
as  temperature  fluctuations.  This  invalidates  many  of  the  assumptions  made  in 
deriving  the  strong  Reynolds  analogy  (SRA).  However,  the  relationship  between  RMS 
temperature  and  streamwise  velocity  fluctuations,  (4.9a),  agreed  with  the  simulation 
data  reasonably  well  nonetheless.  A  condition  for  the  validity  of  the  RMS  relationship 
in  the  presence  of  significant  total  temperature  fluctuations  was  derived,  (4.12),  and 
this  condition  is  satisfied  by  the  simulation  data.  An  expression  for  the  correlation 
coefficient  R^t"  derived  by  Gaviglio  (4.15)  using  the  RMS  relationship  agrees  very 
well  with  the  simulation  data. 

The  low  value  of  the  correlation  coefficient  found  in  the  simulations  indicates  that 
instantaneous  relationships  between  temperature  and  velocity  fluctuations,  (4.8)  for 
example,  are  invalid.  Experimental  evidence,  however,  suggests  a  much  higher  value 
of  the  correlation  coefficient  than  was  found  in  this  simulation.  It  appears  that  this 
difference  between  experiments  and  the  current  simulation  can  be  due  to  a  difference 
of  about  a  factor  of  2  in  the  magnitude  of  the  total  temperature  fluctuations,  with 
the  experimental  values  being  smaller.  The  reason  for  this  is  not  known. 

The  modified  Reynolds  analogy  of  Huang  et  al  showed  better  agreement  with 
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the  simulation  data  than  Gaviglio’s  modified  Reynolds  analogy  and  the  original 
expression  of  Morkovin.  Using  Huang  et  a/.’s  modified  analogy,  a  relationship  between 
the  rate  of  turbulent  heat  transfer  and  turbulent  momentum  transfer  was  derived  and 
shown  to  agree  with  the  simulation  data.  The  streamwise  momentum  and  temperature 
fluctuations  were  found  to  be  very  highly  correlated  throughout  the  boundary  layer 
with  a  correlation  coefficient  0.88  ^  — R^j •  ^  1.  This  is  in  contrast  to  the  low 
correlation  between  the  velocity  and  thermal  fields  away  from  the  wall  and  also  stands 
in  contrast  to  the  lack  of  correlation  between  streamwise  momentum  and  temperature 
in  the  incompressible  case  (where  velocity  and  momentum  are  proportional). 

The  turbulent  kinetic  energy  budget  was  calculated  and  compared  with  those 
of  Spalart’s  incompressible  simulations.  The  peak  rate  of  production  was  found  to 
be  larger  than  for  the  incompressible  case.  This  is  balanced  by  an  increase  in  the 
magnitude  of  turbulent  transport  and  viscous  transport  when  compared  to  the  incom¬ 
pressible  simulations.  Some  of  this  difference  might  be  attributable  to  the  different 
small-scale  resolution  used  in  these  two  simulations,  with  the  current  simulation  hav¬ 
ing  better  resolution  than  the  simulations  of  Spalart.  Balances  for  the  terms  in  the 
Reynolds  stress  tensor  have  been  computed  and  are  presented  in  Guarini  (1998). 

The  authors  would  like  to  thank  NASA’s  Numerical  Aerodynamic  Simulation 
facility  (NAS)  and  the  Air  Force  Aeronautical  Systems  Center  (ASC)  Major  Shared 
Resource  Center  (MSRC)  for  their  computational  support.  The  simulation  and  code 
development  were  done  on  the  Cray  C-90  supercomputers  at  these  two  facilities.  The 
support  of  a  NASA  Graduate  Student  Researchers  Program  grant  (NGT  2-52209) 
and  AFOSR  grant  F49620-97- 1-0089  is  gratefully  acknowledged  by  the  first  and 
second  authors,  respectively. 


Appendix.  Theoretical  development 

In  this  appendix  we  review  Spalart’s  transformation  and  apply  it  to  the  compressible 
boundary  layer.  This  involves  the  development  of  a  generalized  coordinate  system  in 
which  boundary  layer  growth  is  minimal,  the  definition  of  the  two  scales  involved  in 
the  problem,  the  transformation  of  the  Navier-Stokes  equations  to  the  new  curvilinear 
coordinate  system,  and  the  calculation  of  the  slow-growth  terms.  This  analysis  for 
the  compressible  case  mirrors  that  developed  by  Spalart  (Spalart  &  Leonard  1987; 
Spalart  1988)  for  the  incompressible  boundary  layer.  The  particulars  of  that  analysis 
that  are  directly  relevant  to  the  current  development  are  recalled  briefly  in  §  A.2  and 
§  A.3  to  fix  the  ideas  and  the  nomenclature. 

A.l.  Equations 

The  form  of  the  Navier-Stokes  equations  is  chosen  for  computational  convenience. 
The  energy  equation  is  transformed  so  that  pressure  is  a  state  variable  instead  of 
energy.  The  fluid  variables,  u„  m,-  =  pw,-,  p,  and  a  =  (1/p)  are  non-dimensionalized  by 
a0i  p0a0,  p0a20,  and  (l/p0),  respectively.  Lengths  are  non-dimensionalized  by  S0  and 
times  by  ( &0/a0 ).  Then  the  Navier-Stokes  equations  become 


do 

2dmj 

dt 

8xj  ’ 

dmi 

dt 

d  , 

n 

dp 

dpuj 

dt 

dxj 

. 

dX; 


1  dr  u 
Re  dxj 9 


■  + 


1  Hj 
RePr  dx / 


(Ala) 

(Alb) 


(Ale) 
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The  Reynolds  number  is  Re  =  (, p0a0S0)/p0  and  the  Prandtl  number  is  Pr  =  (pCp)/k. 
The  Fourier  heat  conduction  law,  qj,  is  given  by 


«/ 


(A  2) 


and  the  stress,  Ty,  is 

/  5u,-  duj  2  duk 

Tlj  ^  \  dxj  dxi  3  lj  dxk  / 

The  temperature,  T,  is  non-dimensionalized  by  T0  =  such  that  the  equation 

of  state,  T  =  yoP,  results. 


(A3) 


A.2.  Coordinate  system 

Following  Spalart  (Spalart  &  Leonard  1987;  Spalart  1988),  the  Navier-Stokes  equa¬ 
tions  are  transformed  into  a  coordinate  system  that  is  fitted  to  the  growing  boundary 
layer.  The  new  coordinate  system  is  (£,?/, z),  where 


{  =  x  and  rj  =  rj(x,y).  (A  4) 

Curves  of  t]  =  const,  are  slightly  inclined  to  the  wall  with  a  slope  S{^,rj)  chosen  in 
such  a  way  that  we  fit  the  growth  of  both  the  boundary  layer  and  the  viscous  sublayer. 
This  coordinate  system  is  selected  so  that  when  a  small  section  of  the  boundary  layer 
is  simulated,  the  variation  of  the  mean  fluid  dynamic  variables  along  a  constant-?/ 
curve  is  so  small  that  approximate  homogeneity  will  hold. 

The  Jacobian  of  the  transformation  involves  two  parameters,  S  and  7,  where 


s  =  d-Z 

St 


and  7 


=  dy 

dr\ 


(A  5) 


The  quantity  S  gives  the  slope  of  the  constant-?/  curves,  while  7  is  the  local  stretching 
between  the  y-  and  ?/-coordinates,  and  Sn  =  7  follows  from  (A  5).  In  terms  of  S  and 
T,  the  Jacobian  is  given  by 


d/dx  ) 

‘  1 

-S/f 

0  " 

f  d/8Z 

d/dy  >  = 

0 

i/f 

0 

{  d/dr, 

d/dz  ) 

0 

0 

1 

[  d/dz 

(A  6) 


A. 3.  Multiple-scale  analysis 

Even  in  the  transformed  coordinate  system,  the  mean  variables  evolve  slowly  in  <J.  The 
fluctuations  also  have  a  slow  variation  in  intensity  at  constant  rj .  Thus,  for  example, 
we  approximate 

t/(£,  rj ,  z,  t)  =  Au(£9  rj)up{£,  rj ,  z,  f),  (A  7) 

where  Au  is  a  slowly  varying  amplitude  and  up  is  homogeneous,  so  that  in  the 
simulation  it  can  be  treated  as  periodic  in  The  subscript  on  A  refers  to  the 
fluid  dynamic  variable  with  which  it  is  associated.  The  fluctuations  of  the  other 
state  variables  may  be  similarly  decomposed.  After  introducing  a  slow  variable 
E  =  and  a  fast  variable  £,  and  using  the  techniques  of  multiple-scale  asymptotics, 
decomposition  of  the  velocity  into  a  mean  and  fluctuating  part  yields 

u(f ,  E,  rj,  z,  t)  =  U(E9  q)  +  AU(E,  rj)up(Z,  ?/,  z,  t\  (A  8) 


dU  du„  dAu 
e~dE+Au~dZ+eUp^E- 
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and  the  derivative  in  the  streamwise  direction  becomes 

du 

ft 

Using  (A  7)  this  can  be  rewritten  in  the  compact  form 

Wf  =  eUz  +  +  eu3, 

where  i4  —  AUsuf /Au. 

To  allow  u's  =  AUsu'/Au  to  be  determined  in  the  actual  simulation,  note  that 


AuE  _  Kms)g 

A  Kms 


(A  9) 

(A  10) 
t 

(All) 


The  coordinate  system  is  also  slowly  varying  in  £,  and  hence  y  —  y(£,j/),  so  that 
(A  5)  is  rewritten  as 


-4 


and 


T  = 


r\,z 


dy_ 

drj 


(A  12) 


The  simulations  can  be  regarded  as  being  performed  at  a  fixed  value  3  =  S0  of  the 
slow  variable.  We  are  then  free  to  choose  r\  such  that  y(E0,t])  =  t]  which  implies 
T(Eq,  rj)  =  1.  We  define  S  such  that  S  =  eS. 


A.4.  Modified  Navier-Stokes  equations 

Using  the  definitions  above  and  replacing  derivatives  in  £  with  slow  and  fast  deriva¬ 
tives  gives  the  transformed  Navier-Stokes  equations,  that  contain  several  additional 
terms,  shown  enclosed  in  square  brackets: 


do 


7  8m ; 


Sa 


2dmi 

dt]  J 


+  e[—(o3  +  o's)(Tmi+ a(Us +u'3)],  (A  13a) 


dm\  d  dp  1  8xij 

dt  ~  di  +  Red^j+€ 


4(“,m,)+4  j 


+€[{(xs  +  o's)m\  —  2 (Us  +  u'3)mi  -  (p3  +  p's)], 


dm2  d  dp  1  dx2j 

-+T—+t 


S  8 

(«im2  +  w2mi) 


+e[(ffs  +  ff's)mim2  —  (U3  +  u's)m2  —  (Vs  +  i4)»»i], 


dm2  d  dp  i  dx3j 


+e[(ffS  +  a's)m\m3  —  (Ws  +  w4)mi  —  (U3  +  rf3)m3], 
dp  _  dpuj  duj  y-1,  du i  1  dqj 

Tt  ~  ~Ti 7  ~i7~  l)PWj  +  (_R r)XijWj  +  RtfrWj 


+e 


lsA(p„l)+s„-i)p^J 


(A  13  b) 


(A  13c) 


(A  13d) 


+  e[— (Ps  +  p's)°mi  ~  (Us  +  u’s)yp]. 


(A  13e) 
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Note  that  all  additional  viscous  terms  have  been  neglected,  since  they  are  all  multiplied 
by  (S/Re),  which  is  small.  In  particular,  near  the  wall,  where  the  viscous  terms  are 
large,  the  value  of  S  is  approximately  zero.  The  terms  in  square  brackets  are  the 
corrections  to  the  original  Navier-Stokes  equations  that  account  for  boundary  layer 
growth.  In  each  equation  the  first  set  of  bracketed  terms  results  from  the  coordinate 
transformation  and  the  second  set  results  from  the  multiple-scale  analysis.  These 
equations  will  be  solved  in  a  finite  domain  in  the  fast  variable  £.  Thus,  in  the  solution 
domain,  functions  of  the  slow  variable  3  can  be  taken  as  constant  (functions  of  y) 
and  the  fluctuating  quantities  can  be  taken  as  homogeneous  in  the  fast  variable. 


A.5.  Slow  derivatives  of  mean  quantities 

Before  the  modified  equations  (A  13)  can  be  solved  numerically,  the  slow  derivatives 
must  be  determined  in  terms  of  the  simulation  solution  variables.  For  any  arbi¬ 
trary  slowly  varying  function  /,  fs  —  fz/e.  In  what  follows,  it  will  be  convenient  to 
determine  relations  for  rather  than  fz.  The  slow  derivatives  of  the  mean  thermo¬ 
dynamic  quantities  are  calculated  using  the  Van  Driest  (1955)  temperature-velocity 
relationship  as  given  by  Walz  (Fernholz  &  Finley  1980), 


Ml 


(A  14) 


where  the  recovery  factor,  r,  is  taken  to  be  r  =  0.896.  Equation  (A  14)  was  found  to  be 
valid  a  posteriori  in  the  simulations.  Differentiating  (A  14),  the  temperature  derivative 
is  expressed  in  terms  of  the  mean  velocity: 


=  -r(y-l)Ml 


(A  15) 


Since  the  pressure  gradient,  p$,  is  zero  we  get 


Pt 


P 


(A  16) 


Introducing  means  and  fluctuations  into  the  relationship  op  =  1,  averaging,  and 
neglecting  o'pf  gives  after  differentiation 


o 


(A  17) 


Thus  all  the  slow  derivatives  of  mean  thermodynamic  variables  are  related  directly 
to  the  slow  derivative  of  the  mean  streamwise  velocity. 

For  his  incompressible  simulation  at  the  first  station,  Spalart  used  the  well-known 
scaling  laws  for  the  mean  streamwise  velocity  to  calculate  its  slow  derivatives.  There  is 
no  equivalent  scaling  for  compressible  flow.  However,  the  Van  Driest  transformation 
allows  one  to  define  a  transformed  velocity  which  satisfies  the  incompressible  scalings. 
This  transformation  was  found  to  be  valid  a  posteriori  in  the  simulations  (see  §3.1).  In 
the  definition  of  the  van  Driest  transformed  velocity,  Uc  (3.2),  yfTw/T  is  a  function 
of  ^  and  U.  Differentiating  (3.2)  with  respect  to  £  yields 


dUc(Lv)  (V{M  S  (tw\mATr  ,  {Tw\l/2  8U(Z,ti) 

— — l  T({ t)  du+{T> 


(A  18) 


Now,  the  Van  Driest  temperature-velocity  relationship  (A  14)  implies  that  tempera¬ 
ture  is  a  function  of  U  alone  provided  the  recovery  factor,  r ,  is  independent  of 
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Experiments  have  verified  that  this  assumption  is  valid  (Fernholz  &  Finley  1980). 
The  first  term  in  (A  18)  is  therefore  zero  yielding 


To  obtain  UH  (or  U$)  for  use  in  the  modified  Navier-Stokes  equations  (A  13),  the 
strategy  is  to  develop  a  relationship  between  UCi  and  Uc  using  well-known  model 
profiles  for  the  velocity  Uc.  Then  UCi  can  be  calculated  from  Uc  as  determined  in  the 
simulation.  Note  that  model  profiles  discussed  below  are  used  only  to  evaluate  slow 
derivatives .  The  mean  velocity  profile  is  determined  from  the  simulation  not  from  the 
model  profiles. 

To  develop  an  expression  for  UCi,  a  model  profile  for  Uc  that  is  valid  across  the 
entire  boundary  layer  is  needed.  Since  Uc  is  the  Van  Driest  transformed  velocity,  we 
can  use  a  model  profile  for  the  incompressible  boundary  layer,  in  particular,  we  use 
the  relation  of  Coles  (1956): 


UCb/uz  4-  {U/K)w{y/8\  y^8, 
UCoD/uX)  y  >  8, 


(A  20) 


where  UCb  is  a  basic  law-of-the-wall  profile  and  w(y/<5)  is  a  wake  function.  In  order  to 
distinguish  the  model  profile  from  that  of  the  simulation,  the  model  profile  is  denoted 
by  UCm.  For  the  basic  profile  and  wake  function  we  use  the  relations  of  Reichardt 
(1951)  and  Finley  (Cebeci  &  Bradshaw  1977),  respectively: 

^  =  -  In  (l  +  +  C,  [l  -  ,  (A  21) 

Ut  K  \  vw  )  L  V'll vw/  J 


and 


Several  of  the  constants  in  these  expressions  are  prescribed,  that  is  C\  =  — (1/ k:)  In  (k)+ 
C,  7/i  =  ll,  and  b  =  0.33.  The  constants,  uT,  k ,  and  C  are  calculated  from  the  simulation 
mean  velocity  profile  at  each  time  step.  This  leaves  two  parameters,  8  and  FI ,  which 
are  determined  by  matching  the  properties  of  the  model  profile  to  the  instantaneous 
simulation  mean  velocity  profile.  This  profile  is  a  good  representation  of  the  mean 
velocity  throughout  the  boundary  layer  (see  figure  5).  Reichardt’s  profile  is  used 
instead  of  the  classic  profile  of  Coles  because  it  captures  both  the  linear  sub-layer 
behaviour,  where  17+  =  y+,  and  also  the  logarithmic  behaviour  of  the  mean  velocity 
profile.  The  standard  log-law  becomes  infinite  at  the  wall  which  is  undesirable  in  a 
simulation. 

The  parameters,  8  and  FI  are  set  so  that  the  free-stream  velocity  UC(a  and  the 
transformed  displacement  thickness  of  the  model  profile  UCm  match  those  in  the 
simulation.  Thus  we  have 

Uc„  =  UCm(8),  (A23“) 


These  relations  lead  to  a  nonlinear  system  of  equations  for  8 
1998). 


(A  23  b) 
and  IF  (see  Guarini 
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The  slow  derivative  UC;  can  be  written  as 

±k  =  'hi(Ei)  +  JL(¥i) 

«t  «t  \  Mt  /  \  Ux  /  ' 

Using  the  model  profile  to  evaluate  the  derivative,  (A  24)  can  be  written 
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U^  ( Uc\ 

Wt  \Mt  ) 


dUc  8U, 
+  -T7~0(  + 


8b 


dux 


u  +  du*  n 
+  ~mn(’ 


(A  25) 


leaving  just  uXi,  8%,  and  II $  to  be  determined.  It  is  known  that  II  becomes  independent 
of  Reynolds  number  for  Ree  >  5000  (Cebeci  &  Bradshaw  1977)  and  that  the  variation 
of  II  for  lower  Reynolds  numbers  is  very  small,  thus  the  approximation  11^  =  0  is 
used.  Since  Uc^  =  0,  (A  25)  can  be  evaluated  at  y  =  8  to  obtain  a  relation  between 
uxJux  and  8^/8: 


'hi  =  (»t2^/Uc„Vw){(l  +  KbUrM  1  + 1?}  / 

ut  («?<5/l4oVw){(l  +  Kbux/v„)~l  +  R}  +  1  \b  J’ 


(A  26) 


where 


Q~Sux/(v„rn)  _j_  | 

( buxb  t 

]  Q-Suxb/vw 

m 

^  Vw 

) 

(A  27) 


and  the  assumption  S(y  =  <5)  =  b(/b  has  been  made.  To  find  b(/b,  the  momentum 
integral  equation,  and  the  assumption  bjb  =  0^/8  (constant  shape  factor)  is  used  to 
obtain 


b 


Tw 

opJil' 


(A  28) 


This  closes  the  system  of  equations  for  UcJux. 

Since  the  slowly  varying  amplitude  functions  are  proportional  to  the  RMS  fluctu¬ 
ations,  we  can  calculate  the  slow  derivative  of  the  velocity  fluctuations  by  assuming 
a  similar  scaling  law  as  was  used  by  Spalart  for  his  simulations,  which  yields 


Avi  ux 


(A  29) 


Since  the  simulation  results  show  that  the  scaling  used  in  the  incompressible  case 
is  modified  by  the  mean  density  profile,  (A  29)  could  be  improved  by  including  this 
scaling. 

To  determine  the  metric  S  we  take  the  first  derivative  with  respect  to  £  of  the 
following  expression  given  by  Spalart  (1988): 


^(c4y+)  +  f(y/b) 

yp2+yp 


(A  30) 


where  y2  =  {y\y2)'n,  yiux/v„  =  cu  ys/b  =  c2,  and  p  =  c3/logI0(y3/yi).  Spal art’s 
choices  for  the  constants  c\,  C2,  c3,  and  c4  are  0.5,  0.3,  5.0,  and  0.001,  respectively;  rjs  is 
a  weighted  average  of  wall  units  and  y/8  units.  It  should  be  noted  that  the  quantity 
t]s  is  not  the  same  as  y\  and  does  not  satisfy  the  conditions  y  =  q  and  T  ~  1  at  the 
station  of  the  simulation.  Nonetheless,  r\s  can  be  used  to  calculate  S  since  it  follows 
the  growth  of  the  boundary  layer  and  viscous  sublayer  (see  Spalart  1988). 

All  the  slow  terms  in  (A  13)  can  now  be  determined  from  the  simulated  quantities, 
thus  closing  the  equations. 
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THE  NUMERICAL  DECOMPOSITION  OF  TURBULENT 
FLUCTUATIONS  IN  A  COMPRESSIBLE  BOUNDARY  LAYER. 

STANISLAV  G.  BORODAI  AND  ROBERT  D.  MOSER 

Abstract.  In  many  flows  the  turbulence  is  weakly  compressible  even  at  large 
Mach  number.  For  example,  in  a  compressible  boundary  layer  with  Ma  <  5, 
the  differences  relative  to  an  incompressible  boundary  layer  are  understood  as 
being  caused  by  density  variations  that  accompany  variations  in  temperature 
across  the  layer.  Turbulent  fluctuations  in  a  compressible  boundary  layer  are 
therefore  expected  to  be  dominated  by  the  effects  of  non-constant  tempera¬ 
ture,  and  low  Mach  number  theories  in  which  acoustic  fluctuations  are  not 
dominant  should  be  applicable  to  the  fluctuating  field.  However,  the  analysis 
of  compressible  boundary  layer  DNS  data  reveals  the  presence  of  significant 
acoustic  fluctuations.  To  distinguish  between  acoustic  and  thermal  effects,  a 
numerical  decomposition  procedure  for  compressible  boundary  layer  fluctua¬ 
tions  is  applied  to  determine  the  acoustic  and  nonacoustic  fluctuations.  Except 
for  very  near  the  wall,  where  the  decomposition  procedure  is  not  valid,  it  is 
found  that  the  nonacoustic  fluctuations  are  only  weakly  coupled  to  the  acoustic 
fluctuations,  at  Mach  numbers  as  high  as  6. 


1.  Introduction. 

Much  of  the  effort  in  turbulence  research  has  been  directed  towards  incompress¬ 
ible  turbulence,  due  to  the  simplifications  the  incompressible  assumption  provides 
and  because  of  the  belief  that  turbulence  is  primarily  an  incompressible  phenome¬ 
non.  However,  many  flows  of  technological  interest  are  compressible,  so  it  is  of  some 
importance  to  understand  the  effects  of  compressibility  on  turbulence.  Fortunately, 
it  is  observed  in  many  flows  that  turbulence  is  in  fact  only  weakly  compressible, 
even  at  large  Mach  number. 

In  this  paper  we  consider  the  effects  of  compressibility  on  turbulent  boundary 
layers.  A  manifestation  of  weak  compressibility  in  the  boundary  layer  is  Morkovin’s 
hypothesis  [18],  which  suggests  that  the  dynamics  of  a  compressible  boundary  layer 
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is  close  to  incompressible  as  long  as  the  fluctuating  Mach  number  Ma(  is  small. 
Based  on  this  hypothesis,  compressible  boundary  layers  with  free  stream  Ma  up  to 
5  are  considered  weakly  compressible  and  the  differences  relative  to  an  incompress¬ 
ible  boundary  layer  are  understood  as  being  caused  by  the  mean  density  variation 
that  accompanies  the  mean  temperature  variation  across  the  layer.  Morkovin’s  hy¬ 
pothesis  and  other  weakly  compressible  assumptions,  such  as  the  Van  Driest  trans¬ 
formation  [4]  and  Strong  Reynolds  Analogy  [18]  have  been  widely  used  in  modeling 
turbulence  in  compressible  boundary  layers  [26,  13].  Morkovin’s  hypothesis  breaks 
down  when  the  Mach  number  is  so  large  that  the  temperature  fluctuations  are  sub¬ 
stantial  [26].  However,  the  extent  to  which  true  compressibility  effects  (i.e.  finite 
propagation  speed  of  pressure  signals)  are  playing  a  role  in  this  breakdown  is  not 
clear.  If  acoustic  effects  are  not  important  in  this  breakdown,  then  a  formal  low 
Mach  number  asymptotic  treatment  of  the  turbulence  may  be  applicable  to  much 
higher  Mach  numbers. 

1.1.  Low  Mach  number  asymptotic  analysis.  Recently,  a  number  of  papers 
on  low  Mach  number  asymptotic  descriptions  of  compressible  flows  have  been  pub¬ 
lished.  Rigorous  proofs  of  convergence  to  the  incompressible  flow  equations  in  the 
limit  of  low  Mach  number  under  certain  conditions  were  given  by  Klainerman  and 
Majda  [10,  11]  for  Euler  equations  with  polytropic  equation  of  state  and  by  Kreiss 
et  al.  [12]  for  Navier-Stokes  equations.  However,  the  use  of  a  polytropic  equation 
of  state  restricts  the  applicability  of  the  analysis,  by  eliminating  the  heat  transfer 
effects  important  in  the  dynamics  of  a  compressible  boundary  layer.  When  relaxing 
this  restriction,  it  was  pointed  out  by  Matthaeus  et  al.  [28]  and  more  generally  by 
Bayly  et  al.  [1]  that  a  variety  of  low  Mach  number  asymptotic  limits  are  possible, 
depending  on  the  relative  scaling  of  the  thermal  and  acoustic  fluctuations.  Such 
analyses  yield  low  Mach  number  approximations  to  the  compressible  Navier-Stokes 
equations,  which  may  or  may  not  include  acoustic  effects,  depending  on  the  details 
of  the  scaling.  These  approximations  were  also  extended  to  more  complicated  sit¬ 
uations,  such  as  magnetohydrodynamics  [16].  A  number  of  exact  solutions  for  low 
Mach  number  limits  of  the  Navier-Stokes  equations  were  derived  by  Fedorchenko  [6]. 
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Also,  the  low  Mach  number  limit  of  compressible  flow  without  the  assumption  of 
incompressibility  at  the  lowest  order  velocity  was  considered  by  Muller  [19]. 

In  applying  low  Mach  number  asymptotics  to  a  particular  problem,  it  is  neces¬ 
sary  to  determine  which  of  several  Mach  number  scalings  is  appropriate.  For  the 
adiabatic  wall  boundary  layer  turbulence,  thermal  effects  appear  to  dominate  over 
acoustics  since  temperature  fluctuations  are  generated  due  to  the  presence  of  mean 
temperature  gradients.  This  is  consistent  with  approximations  such  as  Morkovin’s 
hypothesis.  The  heat  flux  dominated  hydrodynamics  of  Zank  &  Matthaeus  [28] 
is  one  asymptotic  system  that  appears  appropriate  for  the  boundary  layer.  The 
determination  of  an  appropriate  low  Mach  number  asymptotic  limit  for  boundary 
layer  turbulence  will  provide  a  basis  for  compressible  turbulence  modeling  in  this 
flow,  and  would  provide  a  more  rigorous  foundation  for  such  models  than  the  ap¬ 
proximations  embodied  in  the  Morkovin  hypothesis  and  Strong  Reynolds  Analogy. 

1.2.  Weakly  compressible  turbulence  models.  The  idea  of  utilizing  the  results 
of  weakly  compressible  analysis  in  turbulence  modeling  is  not  new.  It  was  used  ex¬ 
tensively  by  a  number  of  authors.  For  example,  Erlebacher  et  al  [5]  used  a  weakly 
compressible  asymptotic  approach  to  study  the  influence  of  initial  conditions  on 
the  turbulent  dynamics,  Sarkar  et  al.  [23]  and  later  Ristorcelli  [21]  used  it  to  model 
the  dilatational  terms  in  compressible  Reynolds  stress  model,  and  Rubinstein  and 
Erlebacher  [22],  used  this  idea  to  derive  the  transport  coefficients  in  weakly  com¬ 
pressible  turbulence  and,  as  a  result  of  that,  proposed  a  generalized  eddy  viscosity 
model  of  turbulence.  However,  application  of  low  Mach  number  asymptotics  to  tur¬ 
bulence  modeling  is  fraught  with  uncertainties.  In  particular,  one  does  not  know 
a  priori  what  low  Mach  number  scaling  applies  to  the  case  at  hand.  Further,  in 
some  cases,  the  asymptotic  equations  are  written  in  terms  of  distinct  acoustic  and 
nonacoustic  variables.  It  is  very  difficult  to  check  modeling  assumptions  involving 
these  quantities,  since  they  are  not  physically  realized,  and  so  cannot  be  directly 
measured. 

1.3.  Compressible  flow  decompositions.  When  both  acoustic  and  thermal  ef¬ 
fects  are  present,  the  use  of  weakly  compressible  asymptotics  in  modeling  and  other 
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applications  would  be  greatly  facilitated  if  one  could  distinguish  the  acoustic  flue- 
tuations  from  the  turbulence  in  a  compressible  flow  field.  The  idea  of  representing 
a  compressible  flow  as  a  sum  of  different  components  was  introduced  by  Chu  & 
Kovasznay  [3].  They  identified  vortical,  entropic  and  acoustic  components.  How¬ 
ever,  in  their  analysis  Chu  &  Kovasznay  did  not  pursue  the  question  of  how  the 
flow  variables  can  be  actually  decomposed.  Later,  the  mathematical  procedure  of 
decomposition  of  the  flow  variables  was  developed  for  various  flows.  Examples  in¬ 
clude  decomposition  of  disturbances  about  an  arbitrary  potential  mean  flow  [8]  and 
about  transversely  sheared  mean  flow  [7].  The  review  paper  by  Lele  [13]  gives  more 
details  about  these  decompositions.  The  results  of  these  decompositions  may  pro¬ 
vide  certain  advantages  for  the  analysis  of  compressible  flows,  since  the  equations 
for  the  individual  flow  components  are  less  complicated  [20]. 

To  actually  apply  such  a  decomposition  analysis  to  extract  acoustic  components 
of  a  flow,  one  needs  the  three-dimensional  velocity,  density  and  temperature  (or 
pressure)  fields.  Currently,  such  data  are  only  available  from  direct  numerical 
simulation  (DNS).  Such  simulation  data  was  used  by  Blaisdell  [2]  to  decompose 
turbulence  pressure  to  aid  in  the  study  of  the  pressure  strain  correlation.  But,  in 
this  case  only  homogeneous  shear  turbulence  was  considered. 

1.4.  Decomposition  development.  The  research  reported  here  was  undertaken 
to  facilitate  the  use  of  weakly  compressible  asymptotics  for  the  analysis  and  model¬ 
ing  of  turbulence  in  a  compressible  boundary  layer.  Fields  obtained  from  the  direct 
numerical  simulations  of  Maeder  et  al  [14]  at  Mach  numbers  ranging  from  3  to  6 
were  subjected  to  extensive  analysis  to  evaluate  the  validity  of  the  asymptotics  and 
to  characterize  the  compressibility  effects. 

A  preliminary  analysis  (reported  in  section  3)  of  this  data  shows  that  the  HFDH 
of  Zank  &  Matthaeus  [28]  is  an  inadequate  description  because  acoustic  pressure 
and  dilatation  fluctuations  are  not  negligible.  Based  on  these  observations,  a  more 
appropriate  low  Mach  number  scaling  is  determined,  and  governing  equations  for 
acoustic  and  nonacoustic  parts  of  the  flow  are  derived  (section  4).  These  equations, 
extended  to  the  case  of  weakly  compressible  fluctuations  in  the  flow  with  supersonic 
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mean  (appendix  A),  are  used  in  section  5  to  devise  a  numerical  decomposition  algo¬ 
rithm  to  separate  the  acoustic  and  nonacoustic  fluctuations,  and  this  decomposition 
is  validated.  Finally  in  section  6,  the  results  of  the  decomposition  are  discussed, 
along  with  some  concluding  remarks. 


2.  Preliminaries. 

The  research  discussed  here  is  based  on  analysis  of  the  compressible  Navier- 
Stokes  equations,  and  the  direct  numerical  simulation  (DNS)  of  compressible  flow. 


2.1.  Compressible  Navier-Stokes  equations.  In  nondimensional  form  the  com¬ 
pressible  Navier-Stokes  equations  are: 


dp  dpuj  = 

at  dxi  5 

duj  duj  _  1  dp  2  dpsjj 

dt  dxj  'yMa ?  dxi  Re  dxj 

dp  dp  duk  _  7  —  1  d  f  dT\ 
dt  +  U:*  dxj  ^ dx k  Pr  Re  dxk  \  dxk  ) 


+  2 


(7  —  l)7Ma2 
Re 


o  o 

pSij  Sij , 


P  = 


7-1 

7 


pT. 


(1) 

(2) 

(3) 

(4) 


Where 


o  1  f  dui  duj 

Stj  2  VSx.  +  dxi 


j-  r  9Uk 
3  tj  dxk 


is  the  deviatoric  part  of  the  strain-rate  tensor,  V  =  1/p  is  the  specific  volume, 
7  =  cp/cv  is  the  adiabatic  index  with  specific  heats  cv  and  cv  considered  to  be  con¬ 
stants.  The  Fourier  heat  conduction  law  is  used  in  (3)  to  express  the  heat  flux.  The 
second  coefficient  of  viscosity  is  assumed  to  be  zero,  and  the  ideal  gas  law  (4)  is  the 
equation  of  state.  As  is  appropriate  for  boundary  layer  analysis,  free  stream  values 
of  stream- wise  velocity  (uoo)>  density  (p^)  and  pressure  (poo)>  and  the  boundary 
layer  thickness  (5)  were  chosen  as  reference  quantities  to  nondimensionalize  equa¬ 
tions  (l)-(4).  The  coefficients  of  viscosity  p  and  thermal  conductivity  A  were  scaled 
with  respect  to  their  free  stream  values  and  Aqo  ,  and  the  dimensionless  param¬ 
eters  in  (l)-(4)  are  the  Mach  number  Ma  =  u^/coo  where  =  y/^Poo/poo  is  the 
sound  speed,  the  Reynolds  number  Re  =  (poo^oo^)  /  Poo  and  the  Prandtl  number 
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Pr  —  (AiooCp)/A00.  The  Prandtl  number  is  assumed  to  be  constant.  The  dynamic 
coefficient  of  viscosity  ji  is  dependent  on  temperature  T  only,  according  to  the 
Sutherland  law  (see,  for  example,  White  [27]): 


(T\*  I  +  TWTqq 
W  T/a  +  Tsu/Toq  ’ 


(6) 


where  T$u  is  the  Sutherland  constant,  T ^  is  the  free  stream  temperature  and  the 
coefficient  a  =  [T^p^Cp) /p^  =  7/(7—  1)  is  a  consequence  of  the  way  temperature 
was  nondimensionalized  (T  =  T*  p^Cp/poo,  where  T*  is  dimensional  temperature). 

We  use  the  geophysical  coordinate  system,  so  x  is  the  streamwise,  y  is  the  span- 
wise  and  z  is  the  wall-normal  coordinate.  The  (+)  superscript  denotes  the  variables 
in  wall  units,  and  £99  is  used  as  the  boundary  layer  thickness  5.  Reynolds  mean 
quantities  are  obtained  by  averaging  over  streamwise  and  spanwise  directions.  Pro¬ 
files  of  various  quantities  presented  in  the  figures  are  averages  over  time. 


2.2.  Governing  equations  for  fluctuations.  Turbulent  fluctuations  are  the  to¬ 
pic  of  primary  interest.  The  equations  for  the  fluctuations  can  be  obtained  from 
(l)-(4)  by  subtracting  the  equations  for  the  Reynolds  averaged  mean  in  the  usual 
way.  The  result  is: 


dp'  dfj'u,  dpu ( dfh±  \ ' 

dt  dxi  dxj  \  dx{  ) 


7-1  d  (  dT\ 
Pr  Re  dxk  \  dxk  ) 


+  2 


(7  —  Vj'yMa2 
Re 


p'  =  l~  (pT'  +  p'T  +  (p'T'Y)  . 


(7) 

(8) 

(9) 

(10) 


The  overbar  in  the  equations  (7)-(10)  denotes  Reynolds  average,  prime  denotes 
fluctuations.  It  is  common  to  define  the  mean  and  fluctuations  via  the  Favre  average 
(for  a  given  variable  /?,  the  Favre  average  is  /?  =  p/3/p),  since  it  makes  the  mean 
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equations  simpler.  However,  the  use  of  Favre  average  does  not  significantly  simplify 
the  fluctuating  equations,  and  for  the  current  analysis  the  Reynolds  average  was 
found  to  be  more  convenient.  A  similar  analysis  using  the  Favre  average  is  also 
possible. 

2.3.  Properties  of  DNS  data.  In  our  research,  data  from  the  parallel  boundary 
layer  simulations  of  Maeder  et  al  [14]  with  free  stream  Mach  numbers  of  3,  4.5  and  6 
were  used.  The  grid  and  domain  sizes  for  these  simulations  are  shown  in  the  table  1 
(x,  y  and  z  are  stream- wise,  span- wise  and  wall-normal  directions  respectively),  the 
physical  parameters  are  provided  in  table  2. 


Mo 

Grid  size  (Nx  x  Ny  x  Nz) 

Domain  size  ( Lx  x  Ly  x  Lz) 

3.0 

4.5 

6.0 

192  x  144  x  180 

432  x  192  x  200 

240  x  160  x  220 

543+  x  310+  x  1163+ 

1175+  x  282+  x  470+ 

319+  x  193+  x  404+ 

Table  1.  DNS  grid  and  domain  sizes. 


Ma 

Re9 

Tnj/T^ 

Cf  *  103 

Hu 

UT 

3.0 

3028 

2.50 

2.02 

5.86 

0.0498 

4.5 

3196 

4.38 

1.46 

9.22 

0.0552 

6.0 

2945 

6.98 

1.14 

17.20 

0.0614 

Table  2.  DNS  physical  parameters. 


In  Maeder’s  simulations,  the  parallel  boundary  layer  approximation  is  used, 
which  is  based  on  the  assumption  that  the  boundary  layer  is  growing  slowly  in  the 
stream-wise  direction.  This  allows  one  to  consider  a  boundary  layer  with  parallel 
streamlines,  while  boundary  layer  growth  is  accounted  for  by  extra  terms  intro¬ 
duced  into  governing  equations  through  a  two-scale  asymptotic  treatment.  The 
homogenization  technique  used  by  Maeder  is  similar  to  that  of  Spalart  [25]  who 
also  uses  a  two-scale  asymptotic  analysis  in  the  boundary  layer. 

The  horizontal  domain  size  in  Maeder ’s  simulation  is  smaller  than  that  used  by 
other  authors  (e.g.  Guarini  [9]),  so  some  of  the  flow  details  in  Maeder’s  simulation 
may  be  affected  by  domain  size.  However,  the  Maeder’s  data  is  suitable  for  our 
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purposes  since  we  are  primarily  interested  in  compressibility  effects  which  can  be 
investigated  in  small  spatial  domains. 

In  these  simulations,  the  boundary  conditions  at  the  wall  were  no-slip  and 
isothermal  with  temperature  set  to  the  adiabatic  wall  temperature.  For  the  free- 
stream  boundary,  non-reflecting  boundary  conditions  are  used. 


3.  Applicability  of  the  HFDH  approximation. 


Compressibility  effects  observed  in  turbulent  boundary  layers  are  thought  to 
be  dominated  by  variable  fluid  property  effects  (density  and  viscosity)  [26].  One 
would  therefore  expect  the  turbulence  to  be  well  described  by  a  formulation  in 
which  acoustic  effects  play  only  a  minor  role  in  the  flow  evolution,  being  a  higher 
order  correction  in  the  variable  fluid  properties  governing  equations.  In  partic¬ 
ular,  the  Heat  Flux  Dominated  Hydrodynamics  (HFDH)  description  of  Zank  & 
Matthaeus  [28]  would  be  expected  to  be  valid.  To  determine  if  this  is  in  fact  the 
case,  the  DNS  data  of  Maeder  is  examined,  as  described  in  section  2. 

In  the  HFDH  limit,  the  energy  equation  expressed  in  terms  of  the  pressure 
reduces  to  a  compatibility  condition  relating  velocity  divergence  to  the  divergence 
of  heat  flux: 


duk  _  7~I  d  f  dT  \ 

7 Pdxk  Pr  Redxk  \  dxk  )  ’ 


(11) 


Applying  a  similar  analysis  to  the  turbulent  fluctuating  pressure  equation  (9),  one 
obtains  a  relation  for  the  fluctuating  divergence: 


-  dv!k  _  7  -  1  d  (  dT  V 
^ dxk  Pr  Redxk  \dxk) 


(12) 


This  relationship  was  tested  directly  in  the  DNS  data,  and  the  results  are  shown 
in  figure  1.  Throughout  the  boundary  layer,  there  is  a  significant  error  in  the 
satisfaction  of  this  equation;  on  the  order  of  50%.  An  examination  of  the  remaining 
terms  in  the  pressure  equation  shows  that  there  are  several  other  significant  terms 
with  magnitudes  comparable  to  those  in  (12).  The  following  relation  represents 
equation  (12)  extended  to  include  all  of  the  most  significant  terms  in  the  fluctuating 
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0  0.1  0.2  0.3 


r.m.s. 

Figure  1.  The  r.m.s.  profiles  of  the  terms  from  relation 

dv! 

(12):  the  divergence  term  7 ( - ),  the  heat  flux 

term  pfh£-k  {^)  ( . )  and  the  difference  - 

Pr~Re  afr  ( - “)-  The  z  on  this  figure  is  the  wall-normal 

coordinate,  S  =  £99  is  the  boundary  layer  thickness  ( Ma  —  3). 


pressure  equation: 


dpf  _ dp' 

dt  +Ujdxj 


+ 


+  7P 


duk 

dxk 


7-1  d  (  dT  Y 
Pr  Re  dxk  \  dxk  ) 


(13) 


r.m.s.  profiles  of  these  terms  are  presented  in  figure  2.  The  first  three  terms  in  (13) 
(which  are  the  difference  between  (12)  and  (13))  are  large  and  do  not  cancel  each 
other,  despite  the  fact  that  in  the  context  of  the  HFDH  asymptotics,  these  terms 
are  higher  order. 

A  clue  to  the  nature  of  the  discrepancy  is  provided  in  figure  1.  Note  that  outside 
the  boundary  layer  (z/S  <  1),  there  is  significant  velocity  divergence,  but  no  heat- 
flux  divergence.  In  this  region  the  fluctuating  divergence  must  then  be  attributable 
to  acoustics.  The  magnitude  of  this  acoustic  component  of  fluctuating  divergence 
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Figure  2.  The  r.m.s.  profiles  of  the  most  significant  terms  in  flue- 

d\L 

tuating  pressure  equation:  the  divergence  term  ( - ),  the 

heat  flux  term  p~Re  ( . ),  convective  time  deriv¬ 
ative  r  ( - ),  and  convection  with  fluctuating 

velocity  ( - )  (Ma  =  3). 


is  likely  to  stay  the  same  in  the  boundary  layer,  since  there  is  nothing  preventing 
acoustics  from  propagating  into  the  boundary  layer.  Thus,  there  appears  to  be 
an  acoustic  component  of  the  fluctuating  divergence  with  a  magnitude  of  the  same 
order  as  the  magnitude  of  the  divergence  generated  by  the  heat  flux,  suggesting  that 
the  mismatch  in  the  compatibility  condition  (12)  is  due  to  the  acoustic  fluctuations. 

4.  Low  Mach  number  asymptotic  analysis. 

The  goal  of  this  research  is  to  split  compressible  fluctuations  into  acoustic  and 
nonacoustic  parts.  The  governing  equations  for  these  parts  can  be  obtained  con¬ 
sidering  the  limit  of  small  turbulence  Mach  number.  For  simplicity,  asymptotic 
analysis  presented  in  this  chapter  deals  with  the  case  of  low  Mach  number  flow,  the 
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extension  to  the  case  of  low  turbulence  Mach  number  flow  with  supersonic  mean  is 
given  in  Appendix  A.  The  resulting  governing  equations  will  be  used  in  section  5 
to  perform  the  numerical  decomposition  of  compressible  fluctuations. 


4.1.  Asymptotic  scaling.  To  begin  the  analysis,  the  leading  order  asymptotic 
behavior  of  the  flow  variables  as  the  parameter  e  —  yfyMa  -»  0  should  be  postu¬ 
lated,  since  as  described  in  section  1.1  there  are  several  possible  low  Mach  number 
asymptotic  limits.  Based  on  the  observations  of  the  DNS  data  described  in  sec¬ 
tion  3,  the  flow  variables  are  assumed  to  be  composed  of  incompressible  (subscript 
T),  thermal  (subscript  ‘t’)  and  acoustic  (subscript  ‘a’)  parts,  which  scale  with  e  as 
follows: 


u=m  +  e(ut  +  ua), 

(14) 

P=1  +e2(Pi  +Pa), 

(15) 

T=-~  +  eTt  +  s2Ta, 

(16) 

7-1 

p  =1  +  +  £2Pa- 

(17) 

At  zeroth  order  the  flow  variables  are  incompressible,  with  constant  temperature 
and  density  (equal  to  and  1  respectively  because  of  nondimensionalization 
(see  section  2)),  resulting  in  a  constant  zeroth  order  pressure  (equal  to  1  in  this 
nondimensionalization).  However,  the  incompressible  pressure  p/,  which  appears  in 
the  incompressible  Navier-Stokes  equations,  is  of  order  0(e2).  Consistent  with  the 
observation  (section  3)  that  the  acoustic  and  nonacoustic  divergence  are  of  the  same 
magnitude,  the  thermal  and  acoustic  velocity  parts  are  of  the  same  order.  This  fixes 
the  asymptotic  scaling  of  the  other  thermal  and  acoustic  quantities.  Note  that  the 
acoustic  pressure  is  then  the  same  order  as  the  incompressible  pressure,  although 
the  acoustic  density  and  temperature  are  of  higher  order  than  the  corresponding 
thermal  parts.  The  representation  (14)-(17)  is  essentially  an  extension  of  the  HFDH 
limit  of  Zank  &  Matthaeus  [28]  to  include  the  acoustic  component. 
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Since  the  viscosity  depends  only  on  temperature  and  specific  volume  is  just  the 
inverse  of  density,  they  can  also  be  represented  as: 


/i=l  +£Ht+£2Va,  (18) 

V=l  +  eVt  +  e2Va.  (19) 


In  the  limit  of  small  Ma ,  there  are  two  widely  separated  time  scales:  the  in¬ 
compressible  convective  time  scale  5/uoo  and  acoustic  time  scale  (5^/7/coq.  The 
later  time  scale  is  smaller  than  the  former  by  a  factor  of  e.  So  we  can  proceed  as 
in  multi-scale  analysis  and  introduce  two  time  variables,  the  incompressible  (slow) 
time  tj  =  Uoot^/S  =  t  and  the  acoustic  (fast)  time  ta  =  / (5 y/j)  =  t/e .  The 

derivative  with  respect  to  t  is  then  given  by: 


d_  _  ld_  d__ 
dt  e  dta  +  dt]  * 


(20) 


One  can  also  formulate  a  low  Mach  number  asymptotic  analysis  using  multiple 
length  scales.  However,  multiple  length  scale  analysis  is  inappropriate  in  our  case, 
since  the  results  of  this  analysis  will  be  applied  to  the  DNS  data,  and  this  data 
does  not  support  arbitrarily  long  wavelengths. 

Note  that  in  the  current  analysis  we  are  not  using  a  standard  asymptotic  expan¬ 
sion  approach.  In  such  an  approach  all  the  flow  variables  would  be  represented  as 
power  series  of  €  with  order  one  coefficients,  and  then  the  necessary  conditions  for 
the  coefficients  would  be  derived.  Instead,  a  consistent  asymptotic  behavior  has 
been  postulated  through  the  leading  order  scalings  in  (14)-(19).  Further,  the  coef¬ 
ficients  in  (14)— (19)  (e.g.  ut,  ua ,  etc.)  are  order  one,  but  they  are  not  independent 
of  £,  as  would  be  the  case  in  a  standard  asymptotic  expansion.  This  means  that 
there  has  been  no  higher  order  truncation  of  asymptotic  expansion  of  the  variables 
so  far  at  this  point  in  the  analysis.  The  higher  order  terms  in  the  expansions  just 
got  absorbed  into  the  incompressible,  thermal  and  acoustic  variables. 

The  expansions  (14)— (19)  does  not  completely  define  what  the  incompressible, 
thermal  and  acoustic  components  of  the  flow  are.  It  simply  defines  the  leading  order 
scaling  of  these  components.  So,  other  requirements  must  be  imposed  to  uniquely 
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It  follows  from  (22),  that  the  variation  of  nonacoustic  parts  f3tj  on  the  fast  time 
scale  is  two  orders  higher  than  the  leading  order  of  /?*,/,  so  at  the  lowest  two  orders, 
the  nonacoustic  parts  f}tj  can  be  distinguished  from  the  acoustic  part  /3a  by  the  fact 
that  they  only  vary  on  the  slow  time  scale.  This  definition  becomes  ambiguous  at 
higher  orders  at  which  the  nonacoustic  parts  have  fast  time  scale  variation  by  (22), 
therefore  setting  the  maximum  order  that  can  be  considered  in  the  analysis.  The 
highest  order  retained  in  the  equations  for  a  given  variable  /?  will  be  the  order  of 
the  slow  time  derivative  dpa/dtj  of  the  acoustic  part  of  /?. 

Condition  (1)  specifies  the  distinction  between  the  incompressible  and  the  ther¬ 
mal  variables,  while  conditions  (2)  and  (3)  distinguish  the  acoustic  and  non-acoustic 
(thermal  and  incompressible)  quantities  and  provide  a  way  of  separating  governing 
equations  for  these  quantities  at  the  orders  of  interest. 


4.2.  Derivation  of  acoustic  and  nonacoustic  governing  equations.  Our  goal 
is  now  to  formally  decompose  the  governing  equations  into  separate  equations  for 
the  acoustic  and  nonacoustic  parts,  as  in  Chu  &  Kovasznay  [3].  This  is  accomplished 
by  separating  the  terms  that  involve  only  the  slowly  varying  nonacoustic  quantities, 
from  those  that  involve  acoustic  quantities  and  are  therefore  rapidly  varying.  These 
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terms  are  decoupled  on  the  fast  time  scale  since  the  nonacoustic  quantities  do  not 
vary  on  this  time  scale  by  item  (2)  above.  They  are  decoupled  on  the  slow  time  scale 
since  the  fast  acoustic  terms  balance  between  themselves,  even  on  the  slow  time 
scale  (item  3  above).  The  fast  time  scale  equation  will  be  the  acoustic  governing 
equation,  while  the  slow  time  scale  equation  will  be  the  nonacoustic  governing 
equation. 

The  representation  (14)-(20)  is  substituted  into  (l)-(4),  the  terms  of  equal  order 
in  e  are  collected  up  to  the  order  of  the  slow  derivative  of  the  acoustic  quantity,  and 
the  slow  and  fast  terms  are  separated  to  obtain  acoustic  and  nonacoustic  governing 
equations. 

For  example,  consider  the  expansion  of  the  continuity  equation  (1)  (using  (21) 
and  (22)): 


dpt 

dta 


-f  £ 


(dpt  dpa 

\dti  dta 


+  Uij 


dpt  dll ; 


+ 


d'Xj  dxj  dx 


9Ugj\ 
dxj  ) 


2  fdpa  .  dp,  dpt  dpa  dutj 

£  [^  +  Utjd^j+Uajd^  +  U,jd7j+P,'d^  +  Pl 


du 

dx 


) 


+  0(e3)  =  0.  (23) 


The  fast  terms  in  this  equation  yield  the  acoustic  density  equation: 


dpa  ^  duaj 


,  dpa  duai 

+u,id71+pt~di 


,J,aj  ] 

Xj  ) 


+  0(e3)  =  0, 


k  dta  dxj 

while  the  slow  terms  produce  the  nonacoustic  density  equation: 

fdPt 


(24) 


dtj 


Ulj Wi  +  Sf)  +  "2  JhTj  +Pi^~)+  °{e3)  -  °'  (25) 


Applying  the  same  analysis  to  the  the  momentum  equation  (2),  and  keeping 
explicitly  the  terms  up  to  the  order  of  dua/dtj ,  we  get: 


£  dta 
£ 


dun 

dun 

dtj 

+ uij 

dxj 

dUgi 

dun 

dtj 

~+Uti  dxj 

2 

l  _____ 

dsjij 

dxi 

Re 

dxj 

+ 


dUgj  ^ 


dta  dta 


+  u 


aj 


dun 

dxj 


duti 


+  uijd~x ~+Ulj 


duai  \  _ 

f !h~)  ~ 


(26) 


+ 


e  -Vt%L-V&  +  ±\Vt 


dsnj  dptSfij  dstij  ds{ 


dx 


dxi  Re  l  dxj  dxj 


+ 


dxj 


+ 


cnj 


dXj 


+  0{£2)^ 
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where  snj  —  —  % ^/),  &tij  —  —  Lit),  Rud  Sg^j  — -  Siji^u  —  Lia).  Again, 

splitting  (26)  into  acoustic  and  nonacoustic  equations  as  discussed  above,  one  ob¬ 
tains: 


duQ 


dta 


-he 


dUgj 

dtj 


-h  Ujj 


dUgj 

dxj 


dpa 

dxi 


+  6 


Y  ®Pa  .  2  dsgjj 

1  dxi  Re  dxj 


+  0(e2),  (27) 


for  the  acoustic  momentum  equation  and 


dun  dun 

+  ujj  — — -  +  e 


dti 


dxj 

-Vt^-  +  - 
dxi  Re 


du* 


+  utj 


dun 


du- 


dtj  j 


dxj  U! 3  dx 


dpi 

dxi 


2  dsn 


i  \  _  '-'FI  *  WI13 

I  a  i  r»  rv  ' 

i  / 


Re  dxi 


Dsjij  d^tij 

1  dxj  dxj 


+  0(e2)- 


(28) 


dxj  '  dx, 

\  \J 

for  the  nonacoustic  momentum  equation.  As  may  be  expected,  at  lowest  order  the 
equation  (28)  is  just  the  constant  density  incompressible  Navier-Stokes  equation 
governing  u j . 

Finally,  the  procedure  is  applied  to  the  pressure  equation  (3) 


dxj 


dpi  dpa  du 
dtg  +  dta 

.7—1  d2Tt  | 

Pr  Re  dxkdxk 
ji  (  7~1  ( _d_ 


tj 


+^H( 


dpi  dpa  ,  dpi 

Wi+W,  +Ulid^  +  Ulj 


dpA  _ 

dxj  ) 


(29) 


m  .  A 

]  -f 


, Pr  Re  \dxk  \  dxk  J  dxkdxk  J 

yielding  the  acoustic  pressure  equation: 
dpa  duaj\  2  / dpa  dp a  \  _ 

Wa+1^)+£  \dt^  +  Uljd^)- 


\  d2Ta  \  2(7 -1)0  o  \ 

)  +  dxkdxk)+  Re  SliFiij)+0(e  ). 


e2l . }■  . +0(e3), 


Pr  Re  dxkdxk 


(30) 


and  nonacoustic  pressure  equation: 


dutj  2 

ei-Tr^-  +  e 


dx. 


(dpi  ,  .  dpA  _ 
\dti  Ijdxj) 


7-1  d2Tt 

£— — —  — — ^ - H 


i  d 


Pr  Re  dxkdxk 

(  dTt\  2(7 -l)o  o  \ 

\dxk)  +  Sli3Sli3  )  +  )• 


(31) 


KPr  Redxk  dxk  J  '  Re 
At  leading  order,  the  equation  (31)  is  the  analog  of  the  compatibility  condition  (11) 
for  the  nonacoustic  part  of  governing  equations. 

The  leading  order  terms  of  (24), (27)  and  (30)  constitute  the  equations  of  linear 
acoustics.  Further,  by  combining  equations  (24)  and  (30)  one  obtains: 

£dPa  _  _dpa 


dta 


£7C  +  0(1  >■ 


(32) 
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which  implies  that 

e2pa  =  £27  p„  +  0(e3),  (33) 

which  is  just  the  linearized  isentropic  relation  between  pressure  and  density.  In 
(33),  the  additive  0(s2)  function  J{tj,x)  arising  from  time  integration  of  (32)  was 
neglected  by  our  assumption  that  there  is  no  contribution  from  nonacoustic  terms 
varying  on  the  slow  time  scale  to  the  acoustic  equations  at  the  orders  of  interest. 

Now,  for  the  sake  of  completeness,  consider  the  equation  of  state  (4).  The  two 
lowest  nontrivial  orders  are: 

e2(P,  +Pa )  =  -Pt  +  Tt )  + 

7  7-1 

— — - -( PtTt  +  pa — — ~r  +  Ta)  +  0(s3),  (34) 

7  7-1 

Again,  using  the  same  splitting  procedure,  one  obtains  the  acoustic  equation  of 
state: 

e2Pa  =  e2(Pa  +  — irn)  +  0(e3),  (35) 

7 

and  nonacoustic  equation  of  state: 

e2p,  =  e(pt  +  l2Tt)  +  £.2(7  ~  —ptTt  +  0(e3).  (36) 

7  7 

The  leading  order  part  of  nonacoustic  equation  of  state  (36)  is  just  a  Boussinesq 
relation.  Substituting  for  pa  in  (35)  from  (33): 

£2jPa  =  ^2Ta  +  0(e3),  (37) 

which  is  a  linearized  isentropic  relation  between  pressure  and  temperature.  The  re¬ 
lations  (33)  and  (37)  suggest  that  order  e2  acoustic  density  and  temperature  should 
not  be  neglected  in  the  corresponding  asymptotic  expansions,  if  one  is  interested  in 
identification  of  the  acoustic  effects  in  compressible  heat  flux  dominated  flow. 

It  is  instructive  to  contrast  the  above  development  with  the  HFDH  limit  of  Zank 
&  Matthaeus  [28].  Even  though  the  acoustic  pressure  pa  is  included  in  the  HFDH 
analysis,  pressure  fluctuations  do  not  propagate  at  the  first  nontrivial  order.  The 
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reason  is  that  pa  is  neglected  in  HFDH  in  (23).  If  pa  is  neglected,  then  from  (23) 
it  follows: 

dua 


_d 

dt 


!_  ( 

a  \dXj 


+ 


dx  j 


+  0(e)  —  0, 


(38) 


since  Zank  &  Matthaeus  assumed  that  pt  and  uj  do  not  vary  on  the  fast  time  scale 
and  that  uj  is  divergence  free.  If  one  now  takes  the  time  derivative  of  the  equation 
(29)  noting  that  T*  and  p/  do  not  vary  on  the  acoustic  time  scale  in  HFDH,  at 
leading  order  one  gets  (using  (38)): 

02Pa 


dta‘ 


+  0(e)  =  0 


(39) 


rather  then  the  wave  equation  for  pa .  On  the  other  hand,  in  the  current  analysis,  the 
wave  equation  for  acoustic  pressure  can  be  obtained  by  applying  a  similar  procedure 
to  the  equations  (27)  and  (30). 


4.3.  Consolidated  equations.  The  decomposed  equations  derived  in  subsecti¬ 
on  4.2  can  be  rewritten  in  a  more  useful  form.  In  particular,  we  are  interested 
in  the  nonacoustic  components  of  the  variables,  which  include  both  incompressible 
and  thermal  fluctuations.  The  nonacoustic  parts  of  variables  are  therefore  defined 
as  follows: 


un=uI-^eut7  (40) 

pn  =1  +  e2pi,  (41) 

Tn  =  -^-+eTt,  (42) 

7-1 

pn—\-\-ept]  (43) 


where  subscript  V  denotes  the  nonacoustic  part.  The  nonacoustic  equations  are 
then: 


&Pn  ,  dpn'U'nj  ,  r\( ^ 


dt  +  dxj 


Oil ni  vuijii 

~dT  +  Unj~di 7 


du 


+  0(e2)  =  0, 

dUni  1  dpn 


+ 


2  c)  PuS 


mj 


e2pn  dx{  Re  pn  dxj 


+  0(e 2), 


(44) 

(45) 


7  —  1  d 
dxj  Pr  Re  dxk 


nj 


(  dTn\ 


+  0(e2), 


(46) 
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PnTn  =  -2—  +  0{e2).  (47) 

7-1 

We  would  like  to  stress  here  that  all  the  terms  in  the  equations  (44)  (47)  vary  on 
the  slow  time  scale  only.  In  particular,  the  time  derivatives  in  (44)  and  (45)  are 
slow  derivatives.  Also,  it  should  be  noted  that  there  is  an  0(£2)  contribution  to 
nonacoustic  temperature  and  density,  as  a  result  of  the  linear  response  to  the  0(e2) 
incompressible  pressure.  Such  temperature  and  density  variation  is  commonly  re¬ 
ferred  to  as  pseudosound  [17].  The  incompressible  pressure  is  not  dynamic,  and 
so  pseudosound  temperature  and  density  do  not  play  an  important  role  in  the  dy¬ 
namics  of  the  nonacoustic  part  of  the  flow.  They  were  therefore  neglected  in  the 
system  (44)-(47).  That’s  why  equation  (44)  is  written  as  valid  to  the  order  0(£2). 
Note  that  the  nonacoustic  equations  (44)-  (47)  closely  resemble  the  equations  of 
zero  Mach  number  combustion  (see,  for  example  [15]). 

Now,  for  convenience,  the  notation  for  the  acoustic  quantities  is  changed  to  be 
just  the  total  quantity  minus  the  corresponding  nonacoustic  parts,  i.e.: 


Ug- — U  'U'n  — 

(48) 

Ps  =P  -  Pn  -  e2Pa, 

(49) 

TS=T  -  T„  =  e2Tn, 

(50) 

Ps  ~P  ~~  Pn  —  £  Pa] 

(51) 

where  subscript  ‘s’  denotes  the  acoustic  part  in  the  new  notation.  The  correspond¬ 
ing  set  of  governing  equations  for  this  acoustic  part  is: 


dusi  dusi 

+  un 


1  dps 


dt 

dps 


nj  dXj 

dps 


2  dpnsSij 


dt  '  “nj  dXj  •  dXj 


+  unj  Q_  + 

Ps  =  IPs  +  0(e3), 
Ps  —  Ts  +  0(e3). 


e2pn  dxi  Re  pn  dx,j 
dusj  _  7  - 


+  0(e2), 


Pr  Re 


l_d_  (  dTA 

ledxk  \  ndxk) 


+  0(e3), 


(52) 

(53) 

(54) 

(55) 


The  equations  (52)-(55)  are  the  equations  for  acoustic  waves  in  the  convective 
variable  density  medium  defined  by  the  nonacoustic  field. 
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4.4.  Governing  equations  in  the  case  of  a  flow  with  a  supersonic  mean. 
To  study  the  compressibility  effects  in  supersonic  boundary  layers,  the  acoustic 
and  nonacoustic  equations,  derived  above,  must  be  extended  to  the  case  of  a  flow 
with  a  supersonic  mean.  The  extension  relies  on  the  observation  that  the  turbulence 
Mach  number  can  be  small  even  when  the  mean  Mach  number  is  large.  A  low  Mach 
number  asymptotic  truncation  can  therefore  be  fruitfully  applied  to  the  turbulent 
fluctuations,  even  when  it  cannot  be  accurately  applied  to  the  mean.  The  details 
of  this  extended  analysis  are  given  in  the  Appendix  A.  The  resulting  equations  for 
nonacoustic  fluctuations  are  given  by: 


pnTn=^-p  +  0{e2t).  (59) 

7-1 


For  the  acoustic  fluctuations,  we  obtain: 


These  equations  derived  serve  as  the  basis  for  the  numerical  decomposition  of  com¬ 
pressible  boundary  layer  fluctuations  discussed  in  the  section  5. 


5.  Numerical  decomposition. 

To  assist  in  the  evaluation  and  characterization  of  compressibility  effects  in  tur¬ 
bulent  boundary  layers,  a  numerical  decomposition  of  turbulent  fluctuations  into 
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acoustic  and  nonacoustic  parts  is  developed.  The  algorithm  is  based  on  low  turbu¬ 
lence  Mach  number  analysis  discussed  in  A,  and  it  allows  the  decomposed  quantities 
to  be  determined  from  a  compressible  flow  field.  The  description  of  the  algorithm 
and  its  application  to  the  boundary  layer  DNS  data  of  Macder  et  qL  [14]  with  Ma 
equal  to  3,  4.5  and  6  are  described  in  this  section. 

One  of  the  main  goals  of  this  decomposition  is  also  to  directly  evaluate  the  ability 
of  nonacoustic  equations  (56)-(59)  to  describe  the  dynamics  of  the  nonacoustic 
fluctuating  field  in  the  compressible  turbulent  boundary  layer.  These  equations  are 
autonomous,  i.e.  they  do  not  depend  on  the  acoustic  field,  which  depends  on  the 
details  of  the  acoustic  environment  in  which  the  boundary  layer  exists.  So,  if  these 
equations  apply,  they  are  a  valuable  tool  for  the  numerical  modeling  of  compressible 
boundary  layers. 

It  was  noted  when  the  governing  equations  (56) -(59)  for  nonacoustic  fluctua¬ 
tions  were  derived  in  appendix  A,  that  the  time  derivatives  in  these  equations  are 
slow  time  derivatives.  However,  by  our  assertion  that  the  fast  time  derivative  of 
the  nonacoustic  quantity  is  caused  by  acoustic  convection  and  that  these  terms 
cancel  at  the  orders  of  interest  (see  section  A.l),  these  two  terms  can  be  included 
in  the  corresponding  acoustic  or  nonacoustic  governing  equations  and  the  validity 
of  these  equations  will  not  change.  This  will  eliminate  the  problem  of  acoustic  con¬ 
vection  term  being  unaccounted  for  later,  during  the  verification  of  decomposition 
procedure  in  section  5.4.  If  these  terms  are  included  in  the  acoustic  equations,  as 
is  usually  done  in  acoustics,  they  will  increase  the  error  in  the  acoustic  velocity 
equations  by  increasing  the  vortical  part  of  the  acoustic  velocity  time  derivative 
(see  section  5.5  for  details).  Alternatively,  they  are  included  in  the  nonacoustic 
governing  equations.  This  has  the  effect  of  introducing  a  fast  time  derivative  to  the 
nonacoustic  equations.  In  most  cases  (e.g.  turbulence  model  development)  we  are 
only  interested  in  the  equations  on  the  slow  time  scale,  when  (56)-(59)  are  valid. 
However,  when  verifying  the  decomposition  in  section  5.4,  the  total  time  deriva¬ 
tive  will  be  evaluated,  making  it  important  that  the  fast  component  is  included. 
Resulting  equations  will  be  the  same  as  (56)-(59)  with  acoustic  convection  terms 
added. 
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The  derived  governing  equations  for  acoustic  and  nonacoustic  fluctuations  are 
only  valid  and  decoupled  to  a  finite  order,  because  the  assumptions  made  to  separate 
them  are  invalid  at  higher  orders.  Consequently,  there  is  no  information  about  the 
higher  order  coupling  terms,  except  that  their  magnitudes  are  small  relatively  to  the 
magnitudes  of  the  leading  terms  in  the  equations.  No  attempt  will  be  made  to  assign 
these  higher  order  coupling  terms  to  either  acoustic  or  nonacoustic  equations.  After 
the  decomposition  has  been  performed,  these  terms  can  be  computed  explicitly  and 
their  magnitudes  can  be  compared  to  the  magnitudes  of  the  corresponding  leading 
order  terms  to  validate  the  consistency  of  the  approach,  nd  to  assess  the  strength 
of  the  coupling. 

The  numerical  decomposition  is  accomplished  in  stages,  using  the  governing 
equations  for  acoustic  and  nonacoustic  fluctuations.  The  approach  is  to  construct 
fluctuating  fields  that  will  satisfy  the  acoustic  and  nonacoustic  equations  as  closely 
as  possible.  The  chain  of  calculations  used  for  the  decomposition  is  outlined  below. 


5.1.  Temperature  and  density  decomposition.  The  starting  point  of  the  de¬ 
composition  is  the  determination  of  nonacoustic  parts  of  temperature  and  density 
fluctuations.  Since  the  pseudosound  contribution  to  temperature  and  density  was 
neglected  in  the  equations  (56)-(59),  the  nonacoustic  temperature  and  density  fluc¬ 
tuations  for  which  (56)-(59)  are  valid  are  given  by: 


(64) 

(65) 


where  subscript  ‘ps’  denotes  pseudosound  fluctuations.  Pseudosound  fluctuations 
are  defined  to  satisfy  the  linearized  isentropic  relations: 


PPS  (-y-l)T+0^’ 

Ks  =  j  +  otfy, 


(66) 

(67) 


analogous  to  the  isentropic  relations  (62)- (63)  for  acoustic  fluctuations,  with  the  dif¬ 
ference  being  that  pseudosound  fluctuations  are  just  the  linear  response  to  the  non¬ 
propagating  incompressible  pressure  p\ .  Substituting  for  pseudosound  and  acoustic 
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fluctuations  using  (66)-(67)  and  (62)  (63)  from  (64)  the  final  expressions  for  nona¬ 
coustic  density  and  temperature  fluctuations  are  obtained: 


p'n  =  p’  ~ 


P*  _ 
(7-1  )T' 


(68) 

(69) 


So  defined,  nonacoustic  temperature  and  density  satisfy  the  nonacoustic  equation 
of  state  (59)  to  the  order  required,  by  construction.  The  nonacoustic  viscosity 
fluctuations  may  be  obtained  from  Sutherland  law  (6)  with  Tn  =  T  4-  T!n  as  an 
argument: 


i  +  r^/Too  y 
Tn/a  +  Tsu/T0 o  J 


The  nonacoustic  specific  volume  fluctuations  are  defined  as: 


(70) 


(71) 


where  pn  —  p  +  p'n.  Finally,  the  corresponding  acoustic  parts  of  density,  temper¬ 
ature,  viscosity  and  specific  volume  can  be  obtained  by  subtracting  nonacoustic 
parts  from  the  total  fluctuations.  This  will  include  the  pseudosound  component  of 
the  fluctuations  into  the  acoustic  part.  The  pseudosound  component  can  be  sepa¬ 
rated  from  the  acoustic  part  later,  after  the  decomposition  has  been  performed,  if 
desired. 


5.2.  Velocity  decomposition.  The  leading  order  nonacoustic  pressure  equation 


(58): 


-dunj  7 
7 P - -  -  — 


1  d 


dx 


j  Pr  Re  dxk 


(  dTn  V 


2(7  -  l)e2 
Re 


f  o  o  _ o  o' 

pnsijsij  "b  2 pSijSnij 


(72) 


is  just  a  relation  between  nonacoustic  divergence,  nonacoustic  heat  flux  and  leading 
order  dissipation.  Assuming  that  acoustic  dissipation  is  negligible,  the  expression 
for  nonacoustic  divergence  is: 


dxj 


7-1  d  f  dTV\ 

7 p  Pr  Re  dxk  \  n  dxk  ) 


,  2(7  j  l)^2 
7 p  Re 


(73) 


where  Tn  =  T  +  T^,  pn  =  p  +  pfn1  and  non  acoustic  temperature  and  viscosity 
fluctuations  are  obtained  in  5.1.  This  allows  the  divergence  of  the  nonacoustic 
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fluctuations  to  be  determined.  Since  the  divergence  of  nonacoustic  fluctuating 
velocity  is  known,  the  acoustic  divergence  can  be  obtained  by  difference: 


du'Sk  = 

dxk 


du' 


nk 


(74) 

dxk  dxk  dxk 

Now,  apply  the  curl  to  the  leading  order  part  of  the  acoustic  velocity  equation 
(60),  neglect  mean  gradients  and  note  that  the  terms  varying  on  the  slow  time  scale 
don’t  contribute  to  the  acoustic  equations  at  the  orders  of  interest.  The  result  is: 


(^  +  (u  -V))  (V  x  u's)  =  0, 


(75) 


which  is  the  statement  that  leading  order  acoustic  vorticity  is  preserved  in  a  frame 
moving  with  the  mean  flow.  So,  the  leading  order  acoustic  velocity  is  irrotational  at 
all  times,  given  that  it  was  irrotational  at  time  t  —  0.  An  acoustic  velocity  potential 
ip$  can  therefore  be  introduced  to  obtain  a  Poisson  equation  for  <ps  from  (74): 
du ' 


(76) 


Note  that  in  the  presence  of  strong  mean  gradients  the  acoustic  velocity  will  not 
be  irrotational,  so  (76)  will  not  be  valid  close  to  the  wall  in  the  boundary  layer. 

One  also  needs  to  specify  boundary  conditions  for  tps.  Only  the  boundary  con¬ 
ditions  at  the  wall  and  outer  boundary  are  needed,  since  in  the  simulation  data 
being  decomposed,  periodic  boundary  conditions  are  imposed  in  stream-wise  and 
span- wise  directions.  At  the  wall  and  outer  boundary  (z  =  0  and  zmax ): 


9<ps 

dz 


=  Un 


z=0,zn 


(77) 


2=0, -3  rr 


At  the  outer  boundary,  (77)  reflects  the  fact  that  only  acoustic  velocity  fluctua¬ 
tions  are  expected  outside  the  boundary  layer.  At  the  wall,  (77)  is  just  the  no  flow 
through  condition  for  acoustic  and,  consequently,  for  nonacoustic  velocity  fluctua¬ 
tions.  Since  diffusion  of  the  acoustic  field  is  neglected  in  the  acoustic  equations,  the 
no-slip  condition  cannot  be  imposed  on  the  acoustic  velocity  at  the  wall.  This  im¬ 
plies  that  the  tangential  nonacoustic  velocity  will  also  be  nonzero  at  the  wall.  After 
solving  (76)  with  (77)  for  <ps,  one  may  calculate  acoustic  fluctuating  velocity  from 
its  potential  and  nonacoustic  velocity  fluctuations  by  difference,  thus  completing 
the  decomposition  of  the  velocity. 
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5.3.  Pressure  decomposition.  To  decompose  the  pressure  fluctuations  we  will 
seek  nonacoustic  pressure  that  satisfies  the  divergence  constraint,  as  is  done  in 
incompressible  hydrodynamics,  the  difference  being  that  in  our  case  the  divergence 
is  not  zero,  but  given  by  (73).  Taking  the  divergence  of  the  nonacoustic  velocity 
equation  (57)  one  obtains: 


fl¥„  ^i_dvdj/n_ 

dxidxi  V  dx{  dxi 

_  jMa 2  (  d  du'ni  d_  f  du'ni _ 1_  /  ,  dp  ,  dp'n  \  \  V 

V  \dt  dxi  dx-i  \  J  dxj  7 Ma2  \  ndxi  n  dxi  ) ) ) 


(78) 


The  equation  (78)  is  nonlinear  and  can  be  solved  by  iteration.  At  each  iteration, 
essentially  a  Poisson  equation  for  prn  is  solved  with  nonacoustic  specific  volume  on 
the  right  hand  side  of  (78)  computed  with  p'n  from  previous  iteration,  starting  with 
total  pressure  fluctuations  pr  as  initial  guess  for  nonacoustic  pressure.  To  compute 
the  time  derivative  of  the  nonacoustic  divergence,  which  is  part  of  the  right  hand 
side  of  (78),  we  take  the  time  derivative  of  (73)  noting  that  the  nonacoustic  viscosity 
pn  is  a  function  of  Tn  only  when  computing  the  time  derivative  of  viscosity.  The 
Poisson  equation  (78),  when  solved  with  appropriate  boundary  conditions,  gives 
the  nonacoustic  pressure  fluctuation  field. 

The  most  straightforward  boundary  condition  for  p’n  at  the  wall  is  obtained  in  the 
same  way  that  pressure  boundary  conditions  are  obtained  for  incompressible  hydro¬ 
dynamics;  by  evaluating  the  wall-normal  nonacoustic  fluctuating  velocity  equation 
given  by  (from  (57)): 


at  the  wall.  The  first  term  of  the  left  hand  side  of  (79)  is  zero  at  the  wall  because 
of  the  boundary  condition  (77),  and  the  second  term  is  zero  at  the  wall  since  the 
velocity  is  zero  at  the  wall.  The  boundary  condition  for  nonacoustic  pressure  at 
the  wall  is  given  by: 


dPn 

dz 


Z- 0 


2  =  0 


(80) 
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This  boundary  condition  is  valid  provided  the  wall-normal  momentum  equation 
is  satisfied  at  the  wall.  Analytically,  it  is,  but  in  the  numerical  simulations,  the 
governing  equations  are  not  imposed  on  the  boundaries,  instead,  the  boundary 
conditions  are  imposed.  This  leads  to  an  anomaly,  which  is  discussed  in  detail  in 
subsection  5.5. 

The  outer  boundary  condition  is  easier,  because  p'n  is  exponentially  small  there: 
P'n  =  0.  (81) 

Z^Zjnax 

After  prn  has  been  obtained,  it  can  be  subtracted  from  total  pressure  fluctuations  to 
get  acoustic  pressure  fluctuations.  Pseudosound  part  of  temperature  and  density 
can  now  be  separated  from  the  acoustic  part  using  the  isentropic  relations  discussed 
in  section  5.1  if  desired. 

This  step  completes  the  procedure  used  for  the  decomposition  of  a  turbulent 
fluctuating  field  into  nonacoustic  and  acoustic  part.  It  should  be  noted  that  in  the 
homogeneous  case,  by  imposing  the  vorticity  free  and  isentropy  conditions  described 
above  on  the  acoustic  quantities,  one  obtains  fields  that  do  indeed  propagate  as 
acoustic  waves. 

5.4.  Verification  of  decomposition. 

5.4.1.  Acoustic-nonacoustic  pressure  correlation  check.  Asymptotically,  acoustic 
and  nonacoustic  pressure  fluctuations  are  of  the  same  leading  order.  So,  in  princi¬ 
ple,  they  can  be  strongly  correlated.  To  check  that  this  is  not  the  case,  acoustic- 
nonacoustic  pressure  correlation  coefficient  is  computed  and  presented  in  figure  3. 
The  correlation  between  acoustic  and  nonacoustic  pressure  fluctuations  is  weak,  so 
from  this  point  of  view  the  decomposition  is  producing  meaningful  results. 

5.4.2.  Weakly  compressible  nonacoustic  governing  equations  check .  While  develop¬ 
ing  the  decomposition  procedure,  our  goal  was  to  obtain  nonacoustic  variables 
satisfying  derived  weakly  compressible  governing  equations  (56)-(57)  as  close  as 
possible.  The  last  two  of  these  equations  are  satisfied  by  construction,  so  the  valid¬ 
ity  of  the  remaining  two  can  be  checked.  We  can  compute  the  time  derivative  of  a 
given  nonacoustic  variable  /3n  from  the  nonacoustic  governing  equation.  Also,  this 
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Figure  3.  The  profiles  of  acoustic-nonacoustic  pressure  correla¬ 
tion  coefficient  RP'aP'n  for:  Ma  =  3  (- - •),  Ma  =  4.5  ( . ), 

Ma  =  6  ( - ). 


time  derivative  can  be  determined  by  decomposing  two  fields  separated  in  time  by 
small  At  and  computing  the  time  derivative  as  +  At)  -  f3n(At.))/At.  When 
the  time  derivative  of  nonacoustic  variables  is  computed  this  way,  it  includes  the 
the  effect  of  convection  by  acoustic  velocity  (see  section  5).  Therefore,  this  time 
derivative  will  be  compared  to  the  one  from  the  governing  equation  with  the  acous¬ 
tic  convection  included.  If  the  decomposition  is  valid,  then  the  difference  between 
these  two  estimates  of  the  time  derivative  of  will  be  small. 

The  governing  equations  for  acoustic  and  nonacoustic  fluctuations  are  only  valid 
to  a  certain  order.  For  the  velocity  equation  this  is  the  order  of  coupling  terms, 
at  which  the  assumptions  made  to  separate  the  acoustic  and  nonacoustic  equations 
are  no  longer  valid  (see  section  5).  For  the  density  equation,  this  is  the  order  of  the 
explicitly  neglected  pseudosound  terms.  The  sum  of  the  terms  at  these  orders  have 
been  computed  to  provide  an  estimate  of  the  error  committed  by  neglecting  them. 
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The  magnitude  of  these  terms  can  also  serve  as  a  measure  of  how  well  the  nona¬ 
coustic  equations  can  predict  the  evolution  of  turbulent  flow.  Of  course,  the  higher 
order  terms  are  computed  from  the  variables  obtained  by  the  decomposition,  so  this 
measure  will  also  depend  on  the  validity  of  the  decomposition.  The  performance  of 
the  derived  nonacoustic  equations  can  be  assessed  independently  by  performing  a 
DNS  simulation  of  the  flow  governed  by  these  equations  and  comparing  the  results 
to  the  DNS  of  a  fully  compressible  flow.  This  validation,  however,  is  beyond  the 
scope  of  the  current  paper. 

The  r.m.s.  profiles  of  the  time  derivatives  of  nonacoustic  velocity  ufn  computed 
both  ways,  the  difference  between  them  (error),  the  relative  error  and  the  sum  of 
higher  order  terms  in  this  equation  are  presented  in  figure  4.  The  error  is  of  the  same 
order  as  the  sum  of  higher  order  terms,  i.e.  the  error  caused  by  the  approximations 
made  in  the  decomposition  is  of  the  same  magnitude  as  the  intrinsic  error  present 
in  the  equations  due  to  the  neglect  of  higher  order  terms.  The  relative  error  is 
very  small  except  very  close  to  the  wall  (z+  <  30),  where  it  becomes  as  large  as 
20  percent  at  the  wall.  This  is  expected,  since  the  assumptions  made  during  the 
decomposition  become  invalid  close  to  the  wall. 

The  r.m.s.  profiles  of  the  time  derivatives  for  nonacoustic  density  prn  and  the 
errors  are  presented  in  figure  5.  The  behavior  of  the  error  in  the  nonacoustic 
density  equation  is  completely  analogous  to  that  of  the  nonacoustic  velocity. 

It  can  also  be  relevant  to  compare  the  errors  in  the  time  derivatives  to  the  mean 
convective  time  derivatives  D/Dt  —  d/dt  +  ujd/dxj.  The  mean  convective  time 
derivatives  in  our  case  are  factor  of  15  smaller  then  corresponding  time  deriva¬ 
tives.  The  errors  in  the  nonacoustic  equations,  when  compared  to  mean  convective 
derivatives,  are  still  relatively  small. 

5.4.3.  Acoustic  velocity  governing  equation  check .  Here  the  errors  in  the  acoustic 
velocity  equation  (60)  are  presented.  The  order  of  validity  of  the  acoustic  velocity 
equation  is  the  same  as  that  for  the  nonacoustic  velocity  and  it  is  not  influenced 
directly  by  the  neglect  of  pseudosound.  The  profiles  of  time  derivatives  and  errors 
for  acoustic  equation  are  presented  in  figure  6.  The  relative  error  in  this  equation 
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Figure  4.  The  r.m.s.  profiles  of  magnitudes  of  the  time  deriva¬ 
tives  of  nonacoustic  velocity  fluctuations  u'n  calculated  from:  the 

governing  equation  ( - ),  two  subsequent  in  time  data  fields 

( . ),  difference  between  them  x20  ( - ),  magnitude  of 

higher  order  terms  x20  ( - )  and  relative  error  (o  o  o  o)  for: 

a)  Ma  —  3,  b)  Ma  =  4.5,  c)  Ma  =  6. 

is  about  20-30%  in  the  boundary  layer,  but  it  is  very  close  to  the  sum  of  the  higher 
order  terms  throughout  most  of  the  boundary  layer.  The  error  is  also  approximately 
the  same  as  the  error  in  the  nonacoustic  equation.  Relative  error  in  the  acoustic 
equation  is  so  large  simply  because  the  acoustic  time  derivative  is  factor  10  to  20 
smaller  than  the  nonacoustic  time  derivative.  Through  most  of  the  boundary  layer, 
the  errors  in  both  the  acoustic  and  nonacoustic  equations  appear  to  be  inherent 
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Figure  5.  The  r.m.s.  profiles  of  the  time  derivatives  of  nonacous¬ 
tic  density  fluctuations  pfn  calculated  from:  the  governing  equation 
( - ),  two  subsequent  in  time  data  fields  ( . ),  difference  be¬ 
tween  them  x20  ( - ),  magnitude  of  higher  order  terms  x20 

( - )  and  relative  error  (o  o  o  o)  for:  a)  Ma  —  3,  b)  Ma  —  4.5, 

c)  Ma  =  6. 

to  the  splitting  of  the  equation,  rather  then  shortcomings  of  the  decomposition 
procedure,  since  these  errors  are  nearly  the  same  as  the  neglected  higher  order 
terms. 

5.5.  Sources  of  error  in  the  decomposition.  In  summary,  there  are  two  major 
sources  of  error  in  the  decomposition.  The  first,  and  probably  the  most  important, 


30 


STANISLAV  G.  BORODAI  AND  ROBERT  D.  MOSER 


relative  error,  %  relative  error,  % 

0  20  40  60  BO  100  0  20  40  60  80  100 


Figure  6.  The  r.m.s.  profiles  of  magnitudes  of  the  time  derivatives 
of  acoustic  velocity  fluctuations  u's  calculated  from:  the  governing 
equation  ( - ),  two  subsequent  in  time  data  fields  ( . ),  dif¬ 
ference  between  them  ( - ),  magnitude  of  higher  order  terms 

( - )  and  relative  error  (o  o  o  o)  for:  a)  Ma  =  3,  b)  Ma  =  4.5, 

c)  Ma  =  6. 

is  the  neglect  of  the  mean  flow  gradients  when  asserting  that  acoustic  velocity 
is  irrotational,  and  acoustic  fluctuations  are  isentropic.  Because  of  this,  the  time 
derivative  of  acoustic  velocity  us  obtained  from  two  decomposed  consequent  in  time 
data  fields  will  be  irrotational  as  well.  It  is  clear,  however,  that  the  time  derivative 
of  us  as  determined  from  the  governing  equation  (60)  is  not  irrotational.  Vorticity 
is  generated  by  interaction  of  the  acoustic  fluctuations  with  the  mean  gradients,  and 
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this  is  not  accounted  for  in  the  decomposition.  This  inconsistency  is  responsible 
for  the  elevated  error  in  both  the  acoustic  and  nonacoustic  velocity  equations  near 
the  wall. 

The  second  main  source  of  decomposition  error  is  an  anomaly  in  the  DNS  data 
arising  from  the  imposition  of  boundary  conditions  in  the  simulations.  The  de¬ 
composition  described  above  poses  strict  requirements  on  the  DNS  data  used.  In 
particular,  to  avoid  erroneous  boundary  conditions  for  the  nonacoustic  pressure 
p'n  at  the  wall,  which  is  obtained  from  the  governing  equation  for  the  wall-normal 
component  of  fluctuating  nonacoustic  velocity,  we  require  that  the  flow  variables 
satisfy  the  boundary  conditions  used  in  the  simulations,  as  well  as  the  governing 
equations  at  the  wall. 

However,  because  of  the  nature  of  the  numerical  methods  used  to  generate  the 
data  and  the  way  the  boundary  conditions  were  posed  in  the  simulations,  the  gov¬ 
erning  equations  are  not  satisfied  at  the  wall  in  the  DNS  data.  To  overcome  this 
problem,  temperature,  pressure  and  density  were  adjusted  slightly  to  satisfy  the 
governing  equations  at  the  wall.  The  boundary  values  of  pressure,  temperature 
and  density  were  left  unchanged  by  this  data  correction,  only  the  values  at  the 
point  next  to  the  boundary  were  modified.  A  detailed  description  of  the  data 
adjustment  procedure  is  given  in  appendix  B. 

The  efficacy  of  the  data  adjustment  can  be  evaluated  by  examining  the  time 
derivatives  of  the  nonacoustic  pressure  (figure  7).  The  Ma  —  4.5  case  is  presented 
since  it  has  the  biggest  x-y  computational  domain,  which  will  make  the  impact  of 
errors  in  determining  the  pressure  at  the  wall  the  largest.  Here  the  r.m.s.  of  the  pfn 
time  derivative  is  shown  as  determined  by  decomposing  two  fields  separated  by  a 
small  time  increment.  When  the  anomaly  is  not  removed  by  the  data  adjustment, 
there  is  a  large  difference  between  this  time  evolution  and  that  obtained  from  the  p'n 
governing  equation,  but  the  data  adjustment  eliminates  most  of  this  discrepancy. 
Note  however,  that  a  data  adjustment  of  this  sort  is  undesirable;  first  because 
it  is  artificial  and  will  affect  the  subsequent  evolution  of  the  field,  and  second, 
because  even  after  the  data  adjustment  there  is  no  guarantee  that  the  anomaly  has 


32 


STANISLAV  CL  BORODAI  AND  ROBERT  D.  MOSER 


Figure  7.  The  r.m.s.  profiles  of  time  derivatives  of  nonacoustic 

pressure  p'n  computed  from  governing  equation  ( - )  (provided 

for  reference)  and  from  two  subsequent  in  time  data  fields  before 
( - )  and  after  ( . )  the  data  adjustment  (Ma  =  4.5). 


been  eliminated  completely.  So,  DNS  data  that  doesn’t  suffer  from  this  problem  is 
needed. 

The  nonacoustic  fluctuations  defined  by  the  decomposition  satisfy  the  prescribed 
governing  equations.  Also,  the  magnitude  of  higher  order  terms,  reflecting  the 
intrinsic  level  of  error  introduced  by  use  of  weakly  compressible  equations,  is  small. 
This  supports  the  claim  that  the  nonacoustic  flow  evolution  is  autonomous,  i.e. 
it  does  not  depend  on  the  acoustic  field,  and  that  weakly  compressible  governing 
equations  predict  well  the  evolution  of  turbulent  nonacoustic  field.  The  elevated 
error  near  the  wall  in  the  governing  equation  for  acoustic  velocity  is  primarily  due 
to  the  decomposition  deficiencies  discussed  above,  not  the  higher  order  terms,  or 
the  equation  itself.  Efforts  are  under  way  to  improve  the  decomposition. 
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Figure  8.  The  r.m.s.  profiles  of  the  total  ( - ),  acoustic 

( . )  and  nonacoustic  ( - )  fluctuations  of:  a)  velocity, 

b)  pressure,  c)  temperature,  d)  density,  e)  velocity  divergence.  In 

f)  the  total  ( . . --),  acoustic  ( . -)  and  nonacoustic  ( - ) 

Reynolds  shear  stress  is  shown (Ma  =  6). 
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6.  Discussion  and  conclusions. 

The  numerical  investigation  of  the  DNS  data  of  Maeder  et  al  [14]  revealed  the 
presence  of  substantial  acoustic  fluctuations  when  an  application  of  HFDH  asymp¬ 
totics  derived  by  Zank  &  Matthaeus  [28]  to  the  compressible  turbulent  boundary 
layer  fluctuations  was  attempted.  To  study  the  role  these  acoustic  fluctuations  play 
in  the  evolution  of  the  nonacoustic  field,  the  governing  equations  for  the  acoustic 
and  nonacoustic  parts  of  the  flow  were  derived  in  the  limit  of  small  Mach  number. 
These  equations  were  then  adapted  to  the  case  of  small  turbulence  Mach  number;  to 
serve  as  a  basis  of  the  numerical  decomposition  proposed  for  splitting  the  turbulent 
boundary  layer  fluctuations  into  nonacoustic  and  acoustic  parts.  The  numerical  de¬ 
composition  procedure  is  consistent,  in  the  sense  that  the  resulting  acoustic  fields 
satisfy  the  properties  assumed  for  their  derivation,  to  the  order  expected. 

6.1.  Decomposed  quantities.  The  r.m.s.  profiles  of  the  acoustic  and  nonacoustic 
fluctuations  obtained  by  applying  the  decomposition  procedure  to  the  DNS  data 
are  shown  in  figure  8. 

Recall  the  ordering  of  acoustic  and  nonacoustic  fluctuations  postulated  in  sub¬ 
section  4.1  (equations  (14)-(17)),  in  which  acoustic  fluctuations  are  an  order  smaller 
than  the  nonacoustic  fluctuations,  except  for  the  pressure  and  dilatation.  The  r.m.s. 
profiles  shown  in  figure  8  are  consistent  with  this  ordering  throughout  the  boundary 
layer.  But,  outside  the  boundary  layer  (z/5  >  2),  acoustic  fluctuations  dominate  as 
expected.  It  is  remarkable  that  this  ordering  appears  valid  even  at  Ma  =  6,  which  is 
beyond  the  accepted  limits  of  validity  for  weakly  compressible  approximations  such 
as  Morkovin’s  hypothesis.  The  total  {~puw  ~puw),  nonacoustic  (pnU7lwn  —  puw) 
Reynolds  shear  stress  and  the  difference  between  them  associated  with  acoustics 
are  also  shown  in  figure  8,  part  f).  The  acoustic  contribution  to  Reynolds  shear 
stress  is  very  small,  even  at  Ma  =  6.  So  using  the  nonacoustic  equations  to  predict 
Reynolds  stress  in  a  turbulence  model  (say)  is  a  viable  approach. 

6.2.  Nonacoustic  flow  governing  equations.  The  equations  for  nonacoustic 
fluctuations  (56)- (59)  closely  resemble  the  equations  derived  in  HFDH  limit,  how¬ 
ever  in  this  case  these  equations  govern  only  the  evolution  of  the  nonacoustic  part 
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Figure  9.  The  a)  r.m.s.  profiles  of  the  sum  of  higher  order  terms 

( - )  and  of  the  most  significant  of  the  higher  order  terms  in 

density  equation:  u'nj^  ( . ■),  u'sj^  ( - •),  p's ^  ( - ) 

P'slht  (  )  (Ma.  =  6). 

of  the  flow.  Of  these  equations,  the  compatibility  condition  (58)  and  the  equation 
of  state  (59)  are  satisfied  by  construction  of  decomposition  procedure  (see  subsec¬ 
tions  5.2  and  5.1),  while  the  validity  of  the  velocity  and  density  equations  (56)  and 
(57)  was  checked  (see  5.4.2).  The  results  indicate  that  these  equations  predict  the 
evolution  of  the  nonacoustic  field  very  well,  except  near  the  wall,  where  the  error 
appears  to  be  caused  by  the  deficiencies  of  the  decomposition  procedure  for  strong 
inhomogeneities.  The  validity  of  these  equations  at  high  Ma  is  of  great  interest 
since  they  can  be  used  in  turbulence  model  development. 

6.3.  Higher  order  terms.  The  sums  of  higher  order  terms  in  the  fluctuating 
velocity  and  density  equations  were  presented  in  figures  4,  5  and  6.  They  are  also 
shown  in  figures  9  and  10  for  the  highest  Ma  data  studied  (Ma  =  6),  along  with 
the  largest  contributing  terms.  There  are  three  classes  of  higher  order  terms  in 
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Figure  10.  The  a)  r.m.s.  profiles  of  the  sum  of  higher  or¬ 
der  terms  ( - )  and  of  the  most  significant  of  the  higher  or¬ 
der  terms  in  velocity  equation:  ( . ),  ( - ), 

(/inS«u)  (  )  anc^  jfeK  aSj  (^nSn*i)  (  ) 

(Ma  =  6). 


the  fluctuating  equations,  which  were  not  included  in  either  the  acoustic  or  the 
nonacoustic  governing  equations.  They  are: 

1.  The  neglected  terms:  The  convection  of  acoustic  density  with  nonacoustic 
velocity  u'nj  f^j  in  the  density  equation  and  the  viscous  diffusion  of  acous¬ 
tic  velocity  (vnSsij^  *n  fhe  velocity  equation  belong  to  this  class. 

Asymptotically,  these  terms  are  the  most  significant  of  the  higher  order  terms, 
so  they  should  dominate  the  sum  of  higher  order  terms.  This  is  the  case 
for  density  equation.  The  viscous  diffusion  term,  however  is  smaller  than 
other  higher  order  terms  in  the  velocity  equation,  which  is  consistent  with 
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the  assumption  of  acoustic  diffusion  being  negligible.  Also,  since  the  govern¬ 
ing  equations  are  decoupled  at  this  order,  these  terms  clearly  belong  to  the 
acoustic  equations  and  do  not  directly  effect  the  nonacoustic  fluctuations. 

2.  Nonlinear  acoustic  terms:  There  are  relatively  significant  higher  order 
nonlinear  acoustic  terms  u’sj  ^  and  p'$  in  the  density  equation  and  V*8 

in  the  velocity  equation.  Asymptotically,  these  are  the  same  order  as  the 
coupling  terms  discussed  below  and  their  actual  magnitudes  are  smaller  than 
the  dominant  terms. 

3.  Coupling  terms:  These  are  the  higher  order  terms  involving  acoustic  and 
nonacoustic  variables  that  are  present  at  the  order  where  the  governing  equa¬ 
tions  are  no  longer  decoupled.  As  mentioned  in  subsection  5,  there  is  no 
information  about  these  higher  order  terms,  except  for  the  fact  that  their 
magnitudes  are  small  relatively  to  the  magnitudes  of  the  leading  terms  in 
the  equations.  They  cannot  be  assigned  to  either  nonacoustic  or  acoustic 
equations.  These  are  the  true  coupling  terms  that  couple  the  acoustic  and 
nonacoustic  equations.  In  the  density  equation,  the  most  significant  coupling 
term  is  the  interaction  of  acoustic  density  with  nonacoustic  divergence  pls  , 
which  is  of  the  same  magnitude  as  the  nonlinear  acoustic  terms  discussed  in 
item  2,  and  smaller  than  the  dominant  term. 

In  the  velocity  equation,  there  are  two  coupling  terms  that  are  most  sig- 
nificant:  the  nonacoustic  pressure  gradient  term  V8  -gj*  and  the  viscous  diffu¬ 
sion  term  both  caused  by  acoustic  density  variation.  The 

nonacoustic  pressure  gradient  term  is  asymptotically  of  the  same  order  as 
the  nonlinear  term  involving  acoustic  pressure,  discussed  in  item  2.  Its  ac¬ 
tual  magnitude  is  larger  than  one  of  the  nonlinear  term.  This  is  caused  by 
the  difference  in  the  magnitude  of  acoustic  and  nonacoustic  pressure  in  the 
boundary  layer.  The  nonacoustic  pressure  gradient  term  is  dominating  the 
sum  of  higher  order  terms  in  the  velocity  equation,  and  its  relative  magnitude 
is  about  20-30%  of  acoustic  velocity  time  derivative  (see  figure  6). 
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In  judging  the  importance  of  the  neglected  terms  plotted  in  figures  9  and  10,  it 
is  important  to  recall  the  magnitudes  of  the  lowest  order  terms  in  the  equations 
as  shown  in  figures  4  and  5.  Particularly,  the  true  coupling  terms  representing 
interaction  of  acoustics  with  nonacoustic  field  as  discussed  in  item  3  above  are 
quite  small  indeed. 

6.4.  Conclusions.  The  results  presented  in  this  paper  suggest  that  the  nonacous¬ 
tic  equations  (56)-(59)  are  a  valid  description  of  the  evolution  of  nonacoustic  fluctu¬ 
ations  in  a  turbulent  boundary  layer,  up  to  at  least  Ma  =  6.  Since  these  equations 
are  autonomous  (i.e.  do  not  include  acoustic  terms)  they  could  be  used  as  a  sur¬ 
rogate  for  the  compressible  Navier-Stokes  equations,  either  for  model  development 
or  as  a  basis  for  numerical  simulation.  This  would  yield  a  very  good  description  of 
the  turbulent  fluctuations  and  their  effect  on  the  boundary  layer  (e.g.  through  the 
Reynolds  stress).  Another  important  implication  is  that  the  nonacoustic  (i.e.  tur¬ 
bulent)  fluctuations  are  insensitive  to  the  acoustic  fields  present  in  the  flow  under 
consideration.  The  details  of  these  acoustic  fields  depend  on  the  acoustic  environ¬ 
ment  in  which  the  boundary  layer  exists.  Insensitivity  of  turbulence  to  acoustics 
is  necessary  if  boundary  layer  turbulence  at  this  Mach  number  is  to  have  a  mean¬ 
ingful  canonical  state.  Otherwise,  the  structure  of  the  turbulence  will  depend  on 
the  details  of  the  acoustic  environment,  that  is  the  characteristics  of  the  flow  or 
geometry  far  from  the  boundary  layer. 

Our  understanding  of  the  lack  of  coupling  between  acoustic  and  non-acoustic 
fields  in  the  boundary  layer  must  be  qualified  by  the  break-down  of  the  analysis 
very  near  the  wall  where  the  mean  gradients  are  large.  Also,  there  is  a  potential 
coupling  mechanism  right  at  the  wall  since  the  no-slip  boundary  condition  applies 
to  the  sum  of  nonacoustic  and  acoustic  fields,  but  not  necessarily  to  the  fields 
individually. 
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Appendix  A.  Extension  of  low  Mach  number  asymptotic  analysis  to  a 

FLOW  WITH  SUPERSONIC  MEAN. 

In  this  appendix,  the  extension  of  the  low  Mach  number  asymptotic  analysis  of 
section  4  to  the  case  of  the  supersonic  boundary  layer  with  moderately  large  free 
stream  Mach  number  (up  to  say  Ma  =  6)  is  discussed.  Despite  the  supersonic 
free  stream  Mach  number,  if  the  characteristic  turbulence  Mach  number  Mat  = 
Uf/c oo  (where  Uf  is  characteristic  fluctuating  velocity)  of  this  flow  is  relatively  small, 
the  turbulent  fluctuations  can  fruitfully  be  considered  to  be  weakly  compressible. 
In  such  a  case,  an  analysis  analogous  to  that  of  section  4  can  be  applied  to  the 
fluctuating  equations,  with  the  small  parameter  based  on  the  turbulence  Mach 
number. 

There  is  a  subtlety  associated  with  the  application  of  the  asymptotic  analysis  of 
section  4  to  the  fluctuations  in  a  turbulent  boundary  layer.  The  difficulty  occurs  in 
defining  the  zero  Mach  number  limit  being  considered.  The  parameter  describing 
mean  flow  compressibility  is  the  free  stream  Mach  number,  while  the  compressibility 
of  turbulence  is  parametrized  by  the  turbulence  Mach  number.  So,  the  question 
arises  as  to  which  of  the  two  Mach  numbers,  or  both  of  them  should  be  used  as  a 
small  parameter  in  the  asymptotic  analysis  for  weakly  compressible  turbulence. 

In  a  flow  with  nonzero  mean  velocity,  such  as  a  boundary  layer,  the  ratio  of  char¬ 
acteristic  turbulence  Mach  number  to  the  characteristic  mean  flow  Mach  number 
is  equal  to  the  ratio  of  characteristic  velocities,  i.e.  Mat/Ma  =  uj/Uoq .  Therefore, 
considering  the  limit  as  the  turbulence  Mach  number  going  to  zero,  while  the  free 
stream  Mach  number  is  finite  (as  in  Ristorcelli  [21]),  will  require  turbulent  fluctu¬ 
ations  themselves  to  go  to  zero.  However,  the  only  meaningful  limit  is  when  ratio 
Uf/U oo  stays  finite,  requiring  both  free  stream  and  turbulence  Mach  numbers  to  go 
to  zero  together. 

For  the  case  under  consideration,  the  mean  flow  is  supersonic,  while  the  turbu¬ 
lence  Mach  number  is  much  less  the  one.  In  this  case  we  may  expect  an  asymptotic 
truncation  to  be  still  a  good  approximation  for  the  fluctuations,  but  not  for  the 
mean.  Alternatively,  we  could  dispense  with  formal  asymptotics  and  consider  the 
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development  to  be  an  order  of  magnitude  analysis  based  on  the  observed  magni¬ 
tudes  of  fluctuations  in  the  DNS  data. 

One  more  point  worth  discussing  here  is  the  dependence  of  the  temperature  on 
the  small  parameter.  It  is  known  that  for  a  boundary  layer  with  an  adiabatic  wall 
(or  an  isothermal  wall  with  temperature  set  to  the  adiabatic  recovery  temperature 
as  in  the  data  of  Maeder  [14]),  the  mean  temperature  variation  and  the  temperature 
fluctuations  go  to  zero  like  Ma 2  (see  [24]).  This  is  inconsistent  with  the  tempera¬ 
ture  scaling  postulated  for  thermal  fluctuations  in  4.1,  and  with  the  observations 
of  relative  magnitudes  in  the  DNS  data  (see  section  3).  One  way  to  overcome  this 
inconsistency  is  to  consider  a  different  limit  in  which  both  free  stream  and  fluctu¬ 
ating  Mach  numbers  go  to  zero,  but  in  which  the  mean  temperature  and  density 
variations  across  the  boundary  layer  scale  like  Ma  rather  than  Ma2.  For  example, 
the  wall  thermal  boundary  condition  could  vary  with  Ma  to  accomplish  this,  so 
that  the  wall  is  adiabatic  or  at  the  recovery  temperature  only  for  the  Mach  number 
being  analyzed. 

Regardless  of  the  interpretation  of  the  asymptotic  limit,  the  analysis  analogous 
to  that  of  section  4  can  be  applied  to  the  turbulent  fluctuations,  with  small  pa¬ 
rameter  being  the  turbulence  Mach  number.  No  scaling  assumptions  will  be  made 
regarding  the  mean  quantities,  since  the  applicability  of  an  asymptotic  analysis  to 
the  supersonic  mean  flow  is  questionable  for  the  Mach  numbers  of  interest.  Instead, 
the  mean  quantities  will  be  considered  order  one.  To  perform  the  analysis,  the  gov¬ 
erning  equations  (7)-(10)  for  turbulent  fluctuations  need  to  be  rescaled  with  uj  as 
characteristic  velocity. 


A.l.  Scaling.  As  in  section  4.1  the  analysis  is  begun  by  postulating  the  scaling  of 
the  flow  variables  with  the  parameter  et  =  yfyMat: 


U  =  —  U  +  u'l  +  £t{u't  +  u'a), 

(82) 

£/ 

P=P  +  £2i(p'i  +Pa )- 

(83) 

T=T  +  e,T/  +  e\T ", 

(84) 

P=P  +  £tPt  +  £\p'a\ 

(85) 
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The  mean  velocity  scale  is  kept  u0 c,  so  a  factor  e/et  —  Ma/Mat  =  Uoo/uf  appears 
in  front  of  u  in  (82).  Note  that  in  this  analysis  e  =  yfyMa  is  not  a  small  parameter, 
it  is  of  order  one.  The  fluctuations  are  treated  analogous  to  section  4.1,  i.e.  only  the 
leading  order  scaling  is  specified.  The  viscosity  and  specific  volume  can  be  written 
similarly: 


H=H  +  et(j!t  +  £tMo> 

v=v  +  £tv;  +  eX; 

and  the  derivative  with  respect  to  t  is  given  by: 

1A  A 

dt  et  dta  +  dtj ’ 


(86) 

(87) 

(88) 


where  tj  —  Uft*/5  =  t  and  ta  =  Coot* /(by/ 7)  =  t/et. 

Thermal,  acoustic  and  incompressible  parts  of  the  fluctuations  are  defined  to 
have  the  properties  described  in  section  4.1.  The  only  modification  is  that  (22)  is 
modified  to  account  for  the  presence  of  the  mean.  Thus  for  any  quantity  /?,  its 
nonacoustic  parts  fitj  satisfy: 


dAi+£> 

dta  +  1 


,  dW  +  Pji)' 


*aj 


dxj 


u1 


,  dW  +  P'tj)' 


aj 


dxj 


(89) 


Therefore,  fast  time  scale  variation  of  the  nonacoustic  parts  ^[  j  is  primarily  due 
to  convection  of  f3  and  $  7  by  fluctuating  acoustic  velocity. 


A. 2.  Derivation  of  acoustic  and  nonacoustic  governing  equations.  Now, 
by  analogy  with  section  4.2,  we  substitute  (82)-(88)  into  (7)-(10),  collect  the  terms 
of  equal  order  in  et  and  using  the  properties  defined  in  section  A.l  derive  the 
governing  equations  for  acoustic  and  nonacoustic  parts  of  fluctuations.  So,  the 
acoustic  equations  for  density,  velocity  and  pressure  fluctuations,  and  the  acoustic 
equation  of  state  are  then: 
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et 


Re 


£t 


d<i 

dtj 
( 

eV 


4'^y+(^y 


2  I  —dfJ,aSij  ^  cyfdySij  j  ydySgij 


dxj 


a  dxj 


dxj 


+  0(£t)  —  o> 


M  +  eS.M  +  eiv>  *!i  + 

+  £jdxj  +£1Padxj  +lpdx, 


dta 


f) 


*  y<9£/  \  ^dxj)  Pr  Ret  dxk  \  Sx*  adxk) 

— — —  (^VaSijSij  +  2eJiSijSaij\'\  +  0{e3t)  =  0, 


^~Pa  +  ~  ~  ( pT'a  +  Tp'ai'j  +  0{e\)  -  0; 


and  the  nonacoustic  equations  are: 


dp*  ,  duj  ,  dp 

■— +  e^  +  “i  — 


3dx,  ■ 

£'  ((”«fcj)  +  ) +  "t-') = °- 


dn'tj 
+  P~fc~ 


—  eu 


.^4 

j  dxi 


,  dui  ,dp\  da'j  /  ,  9u/vV 

+  £Uljd^  +  Vtfal)  +  ~di7  \Uljlfr~ )  + 


_  ,  <9u*  2 

eUj^7+£Ut^+y^~"jfc; 


cydt-Lt  Sij  +  cyi  dpSjj  +  ydpShj 


dxj 


+  eVl 


•■te+(^),+(^),+«y 


V 


dpstij 

dxj 


+  \Vt 


dxj 


+  V 


Mhii 

dxj 


dxj 
2 

Ret 

+  0(e?)  =  o> 


Sxi 


eV, 


V  Mh 

dxj 


du'tj 

+7^ 


,  dp  (  .  dp  .  duj 

u'^+£t  [u^+nPl^ 

7-1  d  (  dT[  &T\  2(7-1)  ,,2../ 

PrRetdxk\dxk  Ptdxk) 


Re, 


dP'i  .  e 
dti  et 


s  ptSjj$ij  -j-  2epsijSIij  j  j  H“ 

/X  ' 


7-1  »  /  ,dV 


Pr  Ret  dx,k  \ dx^ 
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2^“-  ( ^efiSijStij  +  +  (nhijSiij)  ^  +  °(£t)  =  °> 

et1—  (Tp't  +  pTj)  (-/>'/  +  (P^/)')  +  0(4)  =  0.  (97) 

The  leading  order  terms  of  (90), (91)  and  (92)  constitute  the  equations  for  acous¬ 
tic  fluctuations,  linearized  about  the  mean  flow.  Further,  by  combining  equations 
(90)  and  (92)  one  obtains: 
d  d  duj 


(98) 


£t^£Ujd^+dfa+£dx 


:±)(ek-Pk\  = 

j)  VP  7 Pj 

-  (  Pa  dp  :  p'a  dp  \  ,  2 

£lEUl  V  7 ?dxJ+jPdxj)+°{£t)- 


For  boundary  layer  flow,  the  right  hand  side  of  (98)  can  be  neglected  except  very 
close  to  the  wall,  where  mean  density  gradients  may  be  significant.  So,  from  (98) 
it  follows: 


P'a 


=  +  <*'?>  =  «?(— j)? 


+  0{e\), 


(99) 


which  is  just  the  isentropic  relation  between  pressure  and  density  fluctuations, 
linearized  about  the  mean.  Substituting  for  p'a  in  (93)  from  (99)  one  obtains: 


e\T'a  =  efj:  +  0(4),  (100) 

which  is  the  isentropic  relation  between  pressure  and  temperature  fluctuations, 
linearized  about  the  mean. 


A.3.  Consolidated  equations.  Finally,  the  decomposed  equations  derived  above 
can  be  rewritten  in  a  form  that  will  be  more  useful  in  the  numerical  decompo¬ 
sition  discussed  in  section  5.  They  are  rescaled  with  u0 Q,  and  as  in  section  4.3, 
the  nonacoustic  parts  of  fluctuations  are  introduced  (combined  incompressible  and 
thermal  fluctuations).  For  convenience,  mean  quantities  are  included  in  the  nona¬ 
coustic  variables,  i.e.  nonacoustic  variables  consist  of  the  mean  and  nonacoustic 
parts  of  the  fluctuations.  Also,  as  in  section  4.3,  the  pseudosound  contribution 
to  nonacoustic  temperature,  density  and  divergence  fluctuations  is  neglected.  The 
nonacoustic  equations  are  then  (56)-(59),  presented  in  section  4.4  for  convenience. 
Note  that  the  convective  time  derivative  of  incompressible  pressure  pfj  (terms  in 
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square  brackets  in  (96)),  was  assumed  negligible  in  (58).  This  derivative  is  assumed 
to  be  responsible  for  generation  of  0(e?)  pseudosound  divergence,  analogous  to  the 
way  pseudosound  divergence  is  defined  by  Ristorcelli  [21],  and  should  not  affect  the 
dynamics  of  the  flow  dominated  by  thermal  effects.  Also,  again,  as  in  section  4.3, 
the  time  derivatives  in  (56)--(57)  are  slow  derivatives. 

The  governing  equations  for  acoustic  fluctuations  in  the  new  notation  are  then 
(60)- (63),  presented  along  with  nonacoustic  equations  in  section  4.4.  These  equa¬ 
tions  are  the  equations  for  acoustic  waves  in  the  convective  variable  density  medium 
defined  by  the  nonacoustic  field.  The  acoustic  and  nonacoustic  equations  derived 
in  this  section  will  serve  as  a  basis  for  development  of  the  numerical  decomposition 
procedure  for  boundary  layer  fluctuations  discussed  in  the  section  5. 


Appendix  B.  Data  adjustment  procedure. 


As  discussed  in  section  5.5,  in  order  for  our  decomposition  to  produce  meaningful 
results,  the  DNS  data  must  satisfy  the  governing  equations  at  the  wall  as  well  as  the 
boundary  conditions.  This  is  not  the  case  in  simulations  of  Maeder  et  al  [14]  and 
many  other  DNS,  because  the  governing  equations  at  the  wall  are  usually  discarded 
in  favor  of  imposition  of  the  boundary  conditions.  The  boundary  conditions  used 
in  the  simulation  of  Maeder  et  al.  [14]  are  no-slip  and  isothermal,  so  the  governing 
equations  evaluated  at  the  wall  yield  the  following  relations  for  temperature  and 
pressure: 


d 

(  dT\ 

_ 

dui.  n  O  0  1 

dz 

=  Pr  Re  p- - 27 Ma  PrpSijSjj 

z= 0  ®xk 

dp 

27 Ma2  dfis^j 

dz 

z= 0 

Re  dxj 

2  =  0 

(101) 

(102) 


The  relations  (101)  and  (102)  are  obtained  from  the  thermal  transport  and  the  wall- 
normal  velocity  governing  equations  respectively.  In  this  appendix,  the  procedure 
used  to  adjust  the  DNS  data  to  satisfy  (101)^(102)  is  briefly  discussed. 

In  making  the  adjustment,  the  temperature  at  the  wall  must  remain  unchanged, 
since  it  is  restricted  by  a  boundary  condition.  This  also  means  that  the  viscosity 
at  the  wall  will  remain  unchanged.  To  impose  (101)  then,  the  temperature  at  the 
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first  point  from  the  wall  is  adjusted.  Even  though  there  is  no  boundary  condition 
for  the  pressure  at  the  wall,  the  value  of  pressure  at  the  wall  was  also  chosen  to 
remain  unchanged,  so  the  pressure  is  also  adjusted  at  the  first  point  from  the  wall, 
analogous  to  temperature,  to  satisfy  (102).  Given  the  difference  scheme  used  to 
approximate  the  derivative  at  the  wall,  (101)  and  (102)  become  equations  for  the 
temperature  and  pressure  at  the  first  point  from  the  wall.  These  equations  are 
coupled  due  to  the  presence  of  p,  in  the  derivative  on  the  right  hand  side  of  (102). 
Also,  the  equation  (101)  is  nonlinear  due  to  the  presence  of  p  in  the  derivative  on 
the  left  hand  side.  So,  (101)  is  solved  iteratively  for  the  temperature  at  the  first 
point  from  the  wall,  and  then  pressure  at  the  first  point  from  the  wall  is  determined 
from  (102).  With  new  values  of  pressure  and  temperature,  the  density  at  this  point 
is  determined  from  the  equation  of  state. 

When  this  adjustment  is  applied  to  the  data  of  Maeder  et  al  [14],  the  changes 
made  to  the  temperature  and  pressure  are  relatively  small.  In  the  Ma  =  4.5  case, 
which  is  impacted  the  most  by  the  data  adjustment  procedure,  since  it  has  the 
biggest  computational  domain,  the  r.m.s.  temperature  change  is  5%  of  the  r.m.s. 
temperature  at  the  first  point  from  the  wall,  and  the  r.m.s.  pressure  change  is  less 
then  1%  of  the  r.m.s.  pressure  at  the  first  point  from  the  wall. 

It  should  also  be  mentioned  here  that  the  pressure,  temperature  and  density 
time  derivatives  were  also  adjusted  in  a  similar  manner  using  time  derivatives  of 
conditions  (101)— (102) ,  while  advancing  the  DNS  field  in  time,  to  keep  the  govern¬ 
ing  equations  satisfied  at  the  wall  at  the  next  time  step.  This  was  done  because, 
as  discussed  in  subsection  5.4,  the  decomposition  verification  relies  on  two  decom¬ 
posed  subsequent  in  time  data  fields  to  obtain  the  time  derivatives  of  decomposed 
quantities. 

December  4,  2000 
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Preprint  of  “A  critical  evaluation  of  the  resolution  properties  of  B-spline  and  compact  finite 
difference  methods”  by  W.  Y.  Kwok,  R.  D.  Moser  and  J.  Jimenez  (1999).  Under  consider¬ 
ation  for  publication  in  J.  Comput.  Physics. 
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Resolution  properties  of  B-spline  and  compact  finite  difference  schemes 
are  compared  using  Fourier  analysis  in  periodic  domains,  and  tests  based 
on  solution  of  the  wave  and  heat  equations  in  finite  domains,  with  uni¬ 
form  and  non-uniform  grids.  Results  show  that  compact  finite  difference 
schemes  have  a  higher  convergence  rate  and  in  some  cases  better  resolu¬ 
tion.  However,  B-spline  schemes  have  a  more  straightforward  and  robust 
formulation,  particularly  near  boundaries  on  non-uniform  meshes. 
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1.  INTRODUCTION 

Many  physical  phenomena  involve  a  broad  range  of  spatial  scales.  One  example  is  turbulent 
fluid  flows,  which  have  a  wide  and  continuous  spectrum  of  length  scales  describing  its  compo¬ 
sition  of  eddies  of  different  sizes  [2].  Simulation  of  these  physical  phenomena  requires  spatial 
discretization  schemes  with  high  resolution,  or  in  other  words,  schemes  that  can  produce  accurate 
numerical  results  over  as  broad  a  range  of  length  scales  as  possible  for  a  given  discretization. 

In  numerical  simulation  of  turbulent  fluid  flows,  spectral  methods  are  attractive  spatial  dis¬ 
cretization  schemes  due  to  their  very  good  resolution  properties.  As  a  result,  many  direct  numer¬ 
ical  simulations  (DNS)  have  been  performed  with  spectral  methods  in  Cartesian  coordinates  with 
various  boundary  conditions  (  [4]  and  [14]).  These  include  simulations  of  simple  fundamental 
flows  such  as  isotropic  turbulence,  turbulent  channel  flows  [20]  and  turbulent  boundary  layer  [33]. 
One  distinctive  feature  of  spectral  methods  is  that  they  use  infinitely  differentiable  global  basis 
functions  [4].  Two  common  choices  are  Fourier  series  expansions  and  polynomial  basis  functions, 
with  the  first  being  applied  to  simulations  with  periodic  boundary  conditions  and  the  second  to 
simulations  in  finite  intervals  (  [20]  and  [28]).  However,  the  global  character  of  the  basis  func¬ 
tions  also  limits  spectral  methods  to  simple  geometries  and  boundary  conditions  [23]  and  there 
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is  a  great  need  for  simulations  in  complex  geometries.  This  is  very  important  for  turbulence 
simulations  to  contribute  to  many  engineering  applications  such  as  external  aerodynamics  and 
propulsion  systems.  Such  simulations  would  require  spatial  discretization  schemes  that  not  only 
retain  the  good  resolution  properties  of  spectral  methods,  but  can  also  provide  flexibility  with 
respect  to  geometries,  boundary  conditions  and  grid  distribution. 

Local  numerical  representations,  such  as  finite  difference  and  finite  element  schemes,  have  much 
greater  flexibility  in  discretizing  complex  geometries,  so  high  resolution  schemes  of  these  types 
would  be  of  great  interest.  For  example,  Lele  has  studied  compact  finite  difference  schemes  for 
use  in  problems  with  a  broad  range  of  spatial  scales  [23],  using  Fourier  analysis  to  investigate 
how  well  the  schemes  represent  a  range  of  wave  numbers.  There  has  also  been  a  trend  to  combine 
local  discretization  algorithms  and  spectral  methods.  A  typical  example  of  such  a  confluence  of 
numerical  algorithms  is  the  spectral  element  method  which  is  based  on  finite  element  and  spectral 
methods  [18], [19]  ,[27]. 

Another  choice  for  local  numerical  representation  is  to  use  splines.  Unlike  finite  difference 
methods,  spline  methods  are  functional  expansion  methods  that  make  use  of  a  set  of  local  ba¬ 
sis  functions.  This  property  provides  us  with  a  straight-forward  way  to  implement  boundary 
conditions.  Spline  methods  are  similar  to  finite  element  methods  as  they  both  use  piecewise 
polynomial  representations.  However,  spline  methods  use  basis  functions  that  retain  a  higher 
degree  of  continuity.  In  short,  spline  methods  have  much  of  the  flexibility  afforded  by  the  use 
of  local  expansions,  as  in  finite  elements,  and  have  the  resolution  advantage  afforded  by  highly 
continuous  expansions,  as  in  spectral  methods. 

In  the  research  reported  here,  we  investigate  the  properties  of  spline  methods,  in  particular 
spline  collocation  methods,  and  their  relation  to  finite  difference  and  finite  element  methods.  Sec¬ 
tion  2  introduces  the  basic  properties  of  spline,  compact  finite  difference,  finite  element  methods 
and  their  different  formulation.  The  basic  resolution  properties  of  these  spatial  discretization 
schemes  are  presented  in  Section  3  using  Fourier  analysis  in  periodic  domains.  Of  particular 
interest  is  the  approximation  to  the  first  and  second  derivative  operator,  which  are  common  in 
equations  describing  many  physical  phenomena.  In  Section  4  and  Section  5,  the  wave  equation 
and  heat  equation  are  solved  with  spline  collocation  and  compact  finite  difference  schemes  in 
bounded  domains,  in  both  uniform  and  non-uniform  grids.  Concluding  remarks  are  given  in 
Sections  6. 


2.  NUMERICAL  REPRESENTATIONS 

The  resolution  properties  of  the  numerical  methods  discussed  here  are  most  easily  understood 
in  one  spatial  dimension.  Thus,  the  methods  to  be  evaluated  are  introduced  here  in  their  one¬ 
dimensional  form.  Spline  methods,  compact  finite  difference  methods  and  finite  element  methods 
will  be  discussed. 


2.1.  Spline  Methods 

Consider  a  domain  divided  into  N  intervals,  a  one-dimensional  spline  is  defined  to  be  a  poly¬ 
nomial  of  degree  d  in  each  interval  that  is  continuously  differentiable  d  -  1  times  at  the  interval 
boundaries.  The  boundaries  of  the  intervals  are  called  knots. 

Spline  methods  have  been  used  before  to  solve  differential  equations  and  fluid  mechanics  prob- 
lems(  [13],  [29]  and  [30]).  The  work  of  Kasi  Viswanadham  and  Koneru  [37]  and  Davies  [6], [7] 
used  B-splines  as  basis  functions  and  the  Galerkin  formulation.  Most  of  the  research,  however,  is 
confined  to  cubic  splines  (d  =  3).  More  recently,  Kravchenko  etal.  [21]  and  Shariff  and  Moser  [32] 
used  the  basis  functions  of  splines  to  solve  partial  differential  equations  and  simulate  turbulent 
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FIG.  1.  B-splines  Bf(x)  on  uniform  knots  with  knot  spacing  Ax  =1,  with  (a)  d  =  0,1,2  and  3;  and  with 
(b)  d  =  2  and  i  =  1,2  and  3. 


fluid  flows.  In  particular,  mesh  embedding  techniques  are  developed  to  make  basis  spline  methods 
very  effective  in  solving  physical  problems  in  complex  geometries. 

To  use  splines  as  a  representation  for  the  solution  of  a  partial  differential  equation,  it  is  necessary 
to  have  a  convenient  basis  for  the  space  of  spline  functions  under  consideration.  Here  the  so 
called  basis  splines  or  “B-  splines”  as  described  in  [8]  and  [16]  are  used.  A  B-spline  is  defined  as 
a  normalized  spline  which  has  support  over  the  minimum  possible  number  of  intervals.  In  fact,  it 
has  support  on  only  d  -f  1  intervals.  As  an  example,  the  B-splines  Bf  for  uniformly  spaced  knots 
are  plotted  in  figure  1  (a)  for  d  up  to  3.  By  using  a  basis  with  support  on  the  minimum  possible 
number  of  intervals,  minimum  bandwidth  of  the  resulting  matrices  is  ensured. 

Near  a  boundary,  the  basis  splines  are  different  than  those  in  figure  1(a)  since  the  presence  of 
the  boundary  removes  the  constraint  that  the  B-splines  have  d  —  1  zero  derivatives  at  the  edge 
of  its  interval  of  supports.  An  example  of  the  quadratic  B-splines  near  the  boundary  is  shown  in 
figure  1(6). 

To  use  the  B-splines  in  a  practial  computation,  one  needs  to  evaluate  them  and  their  derivatives 
at  points  in  the  domain.  This  will  be  sufficient  to  compute  the  various  matrices  representing 
different  linear  operators.  An  efficient  and  stable  technique  to  evaluate  the  B-splines  and  their 
derivatives  is  the  recurrence  relation  described  in  [8].  Both  interior  and  boundary  splines  are 
generated  this  way  by  formally  introducing  a  multiplicity  of  knots  at  the  boundary  (see  [8]). 

Consider  the  B-spline  representation  of  a  possibly  non-linear  spatial  operator  T  operating  on 
(f).  We  first  postulate  an  expansion  for  (j)  in  terms  of  B-splines  of  order  don  a  selected  knot  set: 

4>{x)  w  4>{x)  =  ^2aiB?(x)  (1) 


An  approximation  T  to  the  operator  T  is  sought  that  maps  splines  in  Sd  (i.e.  $)  to  splines  in 
Sd ,  where  Sd  is  the  space  of  splines  of  order  d  for  the  selected  knot  set.  That  is 


F(i)Kf(4>)  =  j  =  '£PiB?(x) 

i 


(2) 


There  are  several  ways  to  generate  such  an  approximation.  Two  will  be  considered  here,  namely 
Galerkin  and  collocation  methods. 
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2. LI.  B-spline  Galerkin  Methods 

In  the  Galerkin  formulation,  the  approximation  of  the  linear  differential  operator  V  on  (f>  is 
given  by 

=  (B*,Vj>),  j  =  1,2,...,  (3) 

where  ( f,g )  denotes  the  L2  inner  product  f  fgdx  in  the  domain  and  Nt q  is  the  number  of  B- 
splines.  This  forces  the  error  in  7  to  be  orthogonal  to  Sd,  thus  minimizing  the  L2  error  in  this 
space.  Given  the  linearity  of  V  and  the  representations  of  T>}  4>  and  7,  the  above  equation  can 
be  written 


Nc  Nc 

£  0i(Bf,  Bf)  =  £  atj(Bj  ,T>(Bf))  ,  i  =  l,2 . JVC  (4) 

i= 1  i=  1 

The  inner  products  of  equation  (4)  are  the  elements  of  matrices  M  and  D ,  with  Mij  = 
and  Dij  =  (Bf1V(Bj)).  The  matrix  M  is  called  the  “mass”  matrix  and  D  the  operator  matrix. 
To  obtain  7  given  <f>  one  solves  the  linear  system  M( 3  =  Da.  Note  that  both  M  and  D  are  banded 
matrices  since  individual  B-splines  have  only  local  support.  The  bandwidth  w  of  the  matrices  is 
given  by  w  =  2d  +  1. 

2.1.2.  B-spline  Collocation  methods 

The  collocation  formulation  imposes  different  requirements  to  obtain  the  coefficients  ft.  Here 
the  approximation  7  =  ft#f  of  the  operator  V  on  4>  must  satisfy 

7  =  P^  at  x  =  j  =  (5) 


which  implies 


7VC  Nc 

J^PiBf  =  Y^aiV(Bf)  at  x  =  0,  j  =  1,2, ..  .,N(  (6) 

1=1  1=1 

The  values  of  the  b-splines  and  their  derivatives  are  the  elements  of  the  matrices  M  and  D 
respectively,  with  =  Bf( Q)  and  =  [D(Bf  )](Cj)-  Again,  given  0,  7  is  found  by  solving  the 
linear  system  Mf3  —  Da.  Using  the  collocation  formulation,  the  matrix  bandwidth  w  is  given  by 
w  —  d. 

2.1.3.  Selection  of  Knots  and  Collocation  Points 

To  use  B-splines  in  a  computation,  one  first  needs  to  determine  the  location  of  the  knot  points 
and  for  the  collocation  method  the  collocation  points.  In  a  periodic  domain  with  N  uniform  width 
intervals,  there  are  N  knots  and  N  splines  spanning  the  spline  space.  Therefore,  N  collocation 
points  are  needed  in  a  collocation  scheme.  There  are  only  two  locations  for  the  collocation 
points  that  preserve  the  spatial  symmetry  of  the  operators  :  collocation  points  at  the  knots,  and 
collocation  points  at  the  center  of  the  intervals.  The  former  is  appropriate  for  odd-order  splines, 
the  later  for  even  (see  below). 

In  a  non-periodic  domain,  it  is  more  complicated.  There  are  N  intervals,  N  - j- 1  distinct  knots 
and  N  +  d  collocation  points  are  needed.  There  are  two  basic  ways  to  select  knots  and  collocation 
points  in  a  finite  domain.  The  first  is  that  N  +  d  collocation  points  can  be  selected  by  whatever 
resolution  criteria  are  appropriate  and  then  TV  ~fl  of  these  points  can  be  chosen  to  be  the  knots. 
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Generally  those  collocation  points  that  are  not  knots  are  near  the  boundary,  though  the  knot 
at  the  boundary  is  retained.  This  is  referred  to  as  a  “not-a-knot”  condition,  which  is  commonly 
used  in  spline  interpolation. 

The  alternative  is  to  start  by  selecting  the  knots  according  to  some  resolution  criteria.  This  is 
more  natural  since  the  knots  directly  determine  the  spline  space  and  therefore  are  more  closely 
related  to  resolution  than  the  collocation  points.  Furthermore,  in  a  Galerkin  scheme  all  one  selects 
are  the  knots,  so  direct  comparison  of  Galerkin  and  collocation  is  only  possible  if  one  starts  by 
selecting  the  knots.  Selecting  the  collocation  points  can  then  be  done  in  several  ways,  but  there 
are  two  choices  that  seem  particularly  appropriate  :  place  a  collocation  point  at  the  maximum  of 
each  B-spline  function  or  place  it  at  the  centroid  of  each  B-spline  function.  These  prescriptions 
have  the  advantage  that  they  are  applicable  throughout  the  domain  (nothing  special  about  the 
boundary),  and  they  associate  a  collocation  point  directly  with  each  B-spline  function.  This  later 
property  is  useful  for  applications  in  multidimensional  embedded  grids  of  the  type  described  by 
Shariff  and  Moser  [32].  Note  that  with  uniform  knots  away  from  the  boundary,  the  symmetry 
of  the  B-splines  places  the  maxima  and  centroid  at  the  same  location  :  at  the  knots  or  at  the 
center  of  the  intervals  for  odd  and  even  splines  respectively.  In  the  current  paper,  collocation 
points  at  the  B-spline  maxima  are  selected,  because  this  naturally  places  a  collocation  point  at 
the  boundary,  which  is  useful  for  imposing  boundary  conditions.  Two  knot  distributions  are  used 
:  uniformly  spaced  knots  and  non-uniform  knots  distribution  according  to  : 


x  =  0.5 


cos^Vjf11] } 

«*[*]&■]  j 


(7) 


where  £  =  j /N  for  j  =  0, 1,  •  •  *  N.  This  non-uniform  grid  is  basically  a  Chebyshev  grid  with  the 
boundary  singularities  removed.  It  is  denser  near  the  boundary. 

2.2.  Compact  Finite  Difference  Methods 
Compact  finite  difference  schemes  have  long  been  applied  to  fluid  mechanics  and  other  physics 
problems  [17], [22],  [31].  Recently,  higher  order  compact  finite  difference  schemes  have  seen  in¬ 
creasing  use  in  the  direct  numerical  simulation  of  complex  fluid  flows  (  [12]  and  [26]).  Lele 
presented  a  comprehensive  study  on  the  compact  finite  difference  methods  [23].  Consider  a  uni¬ 
form  mesh  where  the  nodes  are  indexed  by  i.  The  independent  variable  at  the  nodes  is  Xi  and  the 
function  values  at  the  nodes  =  v(xi)  are  given.  The  compact  schemes  are  derived  by  writing 
approximations  of  the  form  : 


pv\_2  +  av[_  i  +  v[  +  av'i+1  +  /?i^+2  —  (8) 

Vi+ 3  —  Vi-3  ,  ,  Vi+2  -  Vi-2  Vi+ 1  —  Vi- 1 

c - — - 1-  b - — - (-  a - — - 

6Arc  4Ax  2Ax 

Similarly,  approximations  to  the  second  derivative  operator  are  derived  by  the  following  rela¬ 
tionship  : 


Pvi- 2  +  avi-i  +  v”  +  av”+  j  +  - 

Ui+3  -  2 Vi  +  Vi-3  ,Vi+ 2  ”  2 Vi  +  Vi-2 

c  9(Aa:)2  4(Az)2 


+  CL 


Ui+1  -  2 Vi  +  Vi-1 


(Ax)2 


(9) 
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TABLE  1 

Table  of  coefficients  for  discretized  first  derivative  operators  using  compact  finite 
difference  schemes.  The  approximations  have  the  form  ^ .  Qj  (fl+j  + 


Band¬ 

width 

Order 

Oi 

OL2 

Ot  3 

04 

Os 

ai 

a2 

a3 

04 

a5 

3 

4 

i 

4 

0 

0 

0 

0 

3 

2 

0 

0 

0 

0 

5 

8 

4 

9 

1 

36 

0 

0 

0 

40 

27 

25 

54 

0 

0 

0 

7 

12 

9 

9 

i 

0 

o 

21 

231 

147 

0 

0 

16 

100 

400 

16 

250 

2000 

9 

16 

16 

4 

16 

1 

0 

144 

152 

10704 

761 

0 

25 

25 

1225 

4900 

125 

125 

42875 

85750 

11 

20 

25 

100 

25 

25 

l 

55 

12760 

5115 

23045 

7381 

36 

441 

784 

15876 

63504 

54 

9261 

10976 

500094 

8001504 

The  relations  between  the  coefficients  a,  6,  c  and  a,/?  are  obtained  by  matching  the  Taylor 
series  coefficients  of  various  orders.  Higher  orders  can  be  obtained  by  including  more  nodes  in 
the  above  two  equations. 

In  this  study,  compact  schemes  with  the  same  stencil  size  on  both  sides  of  the  equations  are 
selected  (c  =  0  in  equations  (8)  and  (9)  for  example).  This  is  because  mass  and  operator  matrices 
with  the  same  bandwidth  is  a  property  shared  by  B-spline  methods.  All  the  coefficients  then 
are  used  to  match  the  Taylor  series  to  as  high  an  order  as  possible.  The  value  of  the  coefficients 
are  listed  in  tables  1  and  2  for  schemes  with  matrix  bandwidth  w  up  to  11.  Note  that  since 
no  restriction  is  imposed  on  the  coefficients  other  than  those  from  Taylor  series  matching,  mass 
matrices  associated  with  the  first,  second  and  higher  derivatives  are  all  different.  This  issue  will 
be  addressed  in  more  detail  in  Section  6. 

2.3.  Finite  Element  Methods 

Most  of  the  finite  element  applications  in  fluid  dynamics  use  the  Galerkin  finite  element  formu¬ 
lation  [11].  The  application  of  the  traditional  finite  element  method  to  fluid  mechanics  is  treated 
by  Thomasset  [35]  and  Baker  [1], 

In  this  study,  the  standard  one-dimensional  Co  finite  elements  are  used.  As  with  B-splines, 
the  finite  elements  are  polynomials  on  a  series  of  knots  (element  boundaries).  However,  because 
only  Co  continuity  is  imposed,  there  are  many  more  degrees  of  freedom  per  interval  (element). 
If  there  are  N  intervals  then  there  would  be  dN  degrees  of  freedom,  where  d  is  the  degree  of 
the  polynomials.  In  this  paper,  only  finite  element  Galerkin  methods  are  considered,  though 
collocation  methods  are  also  possible.  Note  that  this  method  of  increasing  the  local  degree  of 
the  polynomial  shapefunction  is  very  similar  to  the  “p”  finite  element  method  [10],  in  which 
an  element  may  neighbor  an  element  having  different  polynomial  order.  The  main  advantage 
of  finite  element  methods  is  flexibility  with  respect  to  geometry.  Its  weakness  is  resolution. 
In  most  applications  of  finite  element  methods,  elements  are  typically  chosen  to  be  at  most 
quadratic  ([3],  [9]),  and  consequently,  great  accuracy  is  usually  difficult  to  achieve.  This  is 
exactly  opposite  to  the  characteristics  of  spectral  methods.  The  intention  to  combine  these  two 
methods  comprehensively  leads  to  the  development  of  spectral  element  methods  [19],  [27].  Spectral 
element  methods  are  basically  variational  domain  decomposition  techniques.  The  computational 
domain  is  broken  up  into  macro-  elements  within  which  variables  are  represented  as  high-order 
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TABLE  2 

Table  of  coefficients  for  discretized  second  derivative  operators  using  compact 


finite  difference  schemes.  The  approximations  have  the  form  ctj(f”+j  + 


+  Si  — 

Band¬ 

width 

Order 

OL\ 

OL2 

CL% 

OL\ 

C*5 

3 

4 

0.1 

0 

0 

0 

0 

5 

8 

2.91773  X  10”1 

9.75403  x  10“3 

0 

0 

0 

7 

12 

4.54859  X  10_1 

5.07012  x  10~2 

8.53431  X  10"4 

0 

0 

9 

16 

5.70136  x  10”1 

1.12742  x  lO-1 

6.74436  x  10“3 

6.85946  x  10”5 

0 

11 

20 

6.48365  X  10'1 

1.80383  x  10-1 

2.03536  X  10“2 

7.63768  x  10-4 

5.21170  X  10“6 

Band¬ 

width 

Order 

CLl 

O2 

03 

04 

a  5 

3 

4 

1.2 

0 

0 

0 

0 

5 

8 

8.14249  x  lO***1 

7.88804  x  10_1 

0 

0 

0 

7 

12 

3.63508  X  10-1 

1.44465 

2.04670  X  10_1 

0 

0 

9 

16 

7.96376  X  10“2 

1.60567 

6.58960  X  10"1 

3.51163  X  10“2 

0 

11 

20 

-7.02859  X  10“2 

1.47231 

1.12003 

1.72919  X  10_1 

4.77013  X  10~3 

polynomial  expansions  [18].  The  work  of  Patera  [27],  Karniadakis  [18]  and  their  co-workers 
illustrates  the  application  of  spectral  element  methods  in  partial  differential  equations  and  fluid 
mechanics  problems. 


2.4.  Basis  for  Comparison 

To  compare  the  resolution  properties  of  the  several  spatial  discretization  schemes  discussed 
above,  it  is  necessary  to  define  the  basis  of  comparison.  The  question  is  :  comparing  B-spline, 
finite  element  and  finite  difference  methods,  what  characteristics  of  these  methods  (i.e.  what 
degree  polynomials,  or  what  stencil  size)  should  be  compared.  In  this  paper  we  take  the  view 
that  comparison  should  be  done  between  schemes  with  matrices  that  have  the  same  bandwidth. 
The  bandwidth  of  the  matrices  is  an  indicator  of  the  computational  cost  of  the  scheme,  so  methods 
with  the  same  cost  are  compared. 


3.  FOURIER  ANALYSIS 

In  this  section,  a  Fourier  analysis  of  the  errors  associated  with  the  approximation  of  differential 
operators  by  the  several  spatial  discretization  schemes  discussed  in  section  2  is  presented.  The 
resolution  properties  of  the  numerical  schemes  are  most  directly  investigated  using  a  Fourier 
analysis  ([23],  [24],  [25],  [34],  ([36])  ,in  which  the  approximations  of  the  operators  in  a  periodic 
domain  with  a  uniform  grid  are  compared. 
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3.1.  Effective  Wave  Number  and  Eigenfunctions 

One  common  measure  of  how  well  a  differential  operator  is  approximated  is  the  effective 
wavenumber.  In  a  periodic  domain,  the  eigenfunctions  of  derivative  operators  are  the  complex 
exponentials,  and  the  eigenvalues  of  the  nth  derivative  are  (ik)n  where  k  is  the  wavenumber  of  the 
complex  exponential  and  i  ~  y/—l-  The  effective  wavenumbers  k  are  obtained  from  the  eigenval¬ 
ues  of  the  approximate  derivative  operators  M~1D  as  kj  =  ,  where  A  j  is  the  jth  eigenvalue 

of  the  approximate  operators.  For  central  schemes  like  those  studied  in  this  section,  k  is  real. 
Perhaps  more  important  than  the  effective  wavenumber  is  the  error  in  the  eigenvalue  |A  — 

Also  of  interest  is  how  closely  the  eigenfunctions  of  the  approximate  operator  correspond  with 
the  exact  eigenfunctions  (the  complex  exponentials). 

While  the  effective  wave  number  has  been  widely  studied  as  an  indicator  of  the  accuracy 
and  resolution  of  approximate  derivative  operators,  the  accuracy  with  which  the  eigenfunctions 
of  the  operators  (the  complex  exponentials)  are  represented  has  not  generally  been  considered. 
One  reason  is  that  in  finite  difference  methods,  the  circulant  nature  of  the  operator  matrices 
assures  that  the  eigenfunctions  of  the  operators  exactly  represent  the  values  of  etkx  at  the  finite 
difference  grid  points.  However,  with  methods  based  on  functional  representations,  one  can 
measure  the  L2  errors  ||  elkx  —  ^j(x)  ||,  where  rpj(x)  are  the  approximate  eigenfunctions.  Again 
since  the  matrices  M  and  D  are  circulant,  regardless  of  what  derivative  is  being  approximated,  the 
approximate  operator  has  the  same  eigenfunctions  for  all  derivatives.  It  can  be  shown  that  these 
approximate  eigenfunctions  are  also  the  same  as  those  obtained  by  directly  approximating  the 
complex  exponential,  using  the  method  under  considerations  (spline  or  finite  element,  Galerkin 
or  collocation). 


3.2.  Comparison  of  Accuracy  and  Resolution 

The  numerical  schemes  tested  using  Fourier  analysis  include  B-spline  collocation  and  Galerkin 
formulations,  finite  element  Galerkin  formulations,  and  compact  finite  difference  methods.  Effec¬ 
tive  wavenumbers  associated  with  the  first  and  second  derivatives  for  the  four  methods  discussed 
here  are  shown  in  figures  2  and  3  respectively.  Notice  that  the  wavenumber  is  normalized  by  the 
maximum  wavenumber  kmax  representable  with  the  numerical  method.  For  the  spline  and  finite 
difference  methods  kmax  =  3^,  where  Ax  is  the  grid  or  knot  spacing.  For  the  finite  element 
schemes,  kmax  =  since  there  are  d  degrees  of  freedom  per  element. 

There  are  several  things  to  note  about  the  effective  wavenumbers.  First,  for  a  given  matrix 
bandwidth  w,  k  is  identical  for  B-spline  collocation  and  Galerkin  methods.  This  is  despite  the 
fact  that  for  collocation  the  order  of  the  splines  is  higher  ( d  =  w)  than  for  the  Galerkin  (d  =  —p). 
This  identity  had  been  noted  before  by  Swartz  &  Wendroff[34].  Second,  for  a  tridiagonal  matrix, 
the  finite  element  scheme  (linear  elements)  is  identical  to  the  spline  Galerkin  method  (linear 
splines).  For  first  derivatives,  the  effective  wavenumber  is  also  the  same  as  that  for  compact 
finite  difference,  which  is  the  4th  order  Pade  scheme.  For  the  second  derivative,  however,  they 
are  different.  Finally,  the  high-order  (large  bandwidth)  finite  element  effective  wavenumbers 
depart  suddenly  from  the  exact  result,  effectively  limiting  the  range  of  wavenumbers  for  which  k 
is  a  good  approximation  of  k . 

The  errors  in  the  eigenvalues  |A  —  ( ik)n\  for  the  first  and  second  derivatives,  and  those  in  the 
eigenfunctions  are  plotted  in  figures  4,  5  and  6  respectively.  In  comparing  the  different  methods, 
the  most  obvious  difference  is  the  rate  of  convergence  at  small  k :  these  curves  asymptotically 
approach  zero  according  to  their  theoretical  convergence  rate  as  shown  in  table  3. 


fc/fanax  &/fcmax 
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FIG.  2.  Effective  wavenumber  k  of  the  first  derivative  operators  for  matrix  bandwidth  (a)3,  (6)7,  (c)ll: 

- ,  b-spline; - ,  compact  finite  difference; - ,  finite  element  Galerkin; - ,  exact  differentiation. 

For  bandwidth  equals  3,  the  three  schemes  yield  the  same  result. 


k/kmax  k/k  max 
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FIG.  3.  Effective  wavenumber  k  of  the  second  derivative  operators  for  matrix  bandwidth  (a) 3,  (6)7,  (c)  1 1 : 

- ,  b-spline; - ,  compact  finite  difference;  - - ,  finite  element  Galerkin; - ,  exact  differentiation. 

For  bandwidth  equals  3,  b-spline  and  finite  element  yields  the  same  result. 
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fc/^max 

FIG.  4.  Error  in  the  eigenvalue  of  first  derivative  operators  for  matrix  bandwidth  (a)3,  (6)7,  (c)ll:  - , 

b-spline; - ,  compact  finite  difference; - ,  finite  element  Galerkin.  For  bandwidth  equals  3,  all  schemes 

yield  the  same  result. 

Note  that  the  compact  finite  difference  convergence  rate  is  significantly  faster  for  large  w.  This 
is  possible  because  in  the  finite  difference  method,  the  “mass”  matrix  can  be  different  for  each 
order  derivative.  In  contrast,  by  the  nature  of  functional  expansion  methods,  the  mass  matrix  is 
the  same  for  all  derivatives  that  can  be  determined  from  the  representation.  If  this  restriction 
were  imposed  on  the  compact  finite  difference  methods,  the  same  order  of  convergence  as  the 
spline  and  finite  element  methods  would  be  obtained. 


TABLE  3 

Order  of  convergence  of  the  errors  of  eigenvalues  and  eigenfunctions.  w  is 

the  matrix  bandwidth. 


Numerical 

Scheme 

Eigenvalue 
of  the  first 
derivative 
operator 

Eigenvalue 
of  the  second 

derivative 

operator 

Eigenfunction 

finite  element  Galerkin 

kw+ 2 

kw+i 

t =** 

B-spline  Galerkin 

kw+2 

kw+i 

w  +  l 

k 2 

B-spline  collocation 

kw+2 

kw+ 1 

kw+1 

compact  finite  difference 

k2w-l 

k2w 

N.A. 

The  other  property  of  the  approximate  operators  is  the  behavior  of  the  error  at  large  k.  This  is 
important  because  it  determines  the  range  of  spatial  scales  that  can  be  resolved  by  the  numerical 
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FIG.  5.  Error  in  the  eigenvalue  of  second  derivative  operators  for  matrix  bandwidth  (a)3,  (6)7,  (c)ll: 

- ,  b-spline; - ,  compact  finite  difference; - ,  finite  element  Galerkin.  For  bandwidth  equals  3, 

b-spline  and  finite  element  yields  the  same  result. 


method.  There  is  no  commonly  used  measure  of  this  resolution  property  of  numerical  methods. 
One  measure  proposed  by  Lele  [23]  is  the  lowest  wavenumber  (k/kmax)  at  which  the  error  crosses 
some  arbitrary  threshold  (say  0.1),  giving  the  fraction  of  the  maximum  wavenumber  range  that 
is  represented  to  this  accuracy  or  better.  In  table  4,  this  resolved  fraction  for  10  %  error  in  the 
eigenvalues  and  eigenfunctions  is  listed  for  the  numerical  schemes  discussed. 

Despite  the  fact  that  the  order  of  convergence  for  the  finite  element  and  spline  effective  wave 
numbers  is  the  same,  the  errors  in  the  spline  methods  are  much  lower  at  any  given  k.  In  essence, 
the  spline  methods  have  better  resolution.  This  is  indicated  by  Lele’s  resolution  measure  as 
shown  in  table  4.  The  reason  for  the  relatively  poor  resolution  of  the  finite  element  is  the  low 
(Co)  continuity  at  the  element  boundary.  One  way  to  understand  this  (for  the  first  derivative)  is 
to  imagine  a  high  order  finite  element  function  u  evolving  according  to  the  scalar  wave  equation: 
W  +  c§^  ~  0.  At  the  initial  time  there  are  discontinuities  in  first  derivative  at  the  element 
boundaries.  The  exact  solution  would  have  these  discontinuities  propagate  into  the  middle  of  the 
element,  where  they  cannot  be  well  represented,  leading  to  relatively  large  errors.  This  scenario 
suggests  that  maximum  possible  continuity  at  the  knots,  that  is  splines,  is  desirable. 

The  uniform  grid  periodic  analysis  is  informative,  but  it  does  not  address  two  key  issues 
commonly  encountered  in  numerical  simulations,  that  is  nonuniform  grids  and  boundaries.  The 
behavior  of  finite  difference  methods  in  particular  is  at  issue  since  the  formulation  discussed  in 
section  2.2  does  not  apply  directly  in  these  cases.  Also,  the  result  that  the  eigenfunctions  of  the 
derivative  operators  are  recovered  exactly  (section  3.1)  will  not  hold.  It  is  thus  of  interest  to 
consider  model  problems  in  finite  domains.  Two  such  problems  are  discussed  here,  namely  the 
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FIG.  6.  1/2  error  in  the  representation  of  the  complex  exponential  with  wavenumber  k  for  matrix  bandwidth 

(a)3,  (6)7,  (c)ll:  - ,  b-spline  collocation; - ,  b-spline  Galerkin; - ,  finite  element  Galerkin.  For 

bandwidth  equals  3,  b-spline  Galerkin  and  finite  element  Galerkin  yields  the  same  result. 

TABLE  4 

10  %  resolved  fraction  for  eigenvalues  of  first  and  second  derivative  operators, 

and  eigenfunction  representation. 


d/dx  d2/dx2 


Band¬ 

width 

B-spline 

Compact 

finite 

difference 

Finite 

element 

Galerkin 

Band¬ 

width 

B-spline 

Compact 

finite 

difference 

Finite 

element 

Galerkin 

3 

0.59 

0.59 

0.59 

3 

0.68 

7 

0.80 

0.84 

0.82 

7 

1.00 

0.94 

0.57 

11 

0.87 

0.90 

0.54 

11 

1.00 

0.99 

0.59 

Eigenfunction 

Band-  B-spline  B-spline  Finite 

width  collocation  Galerkin  element 


Galerkin 


3 

0.68 

0.46 

0.46 

7 

0.84 

0.72 

0.48 

11 

0.89 

0.81 

0.53 

wave  and  heat  equation.  They  will  only  be  applied  to  the  B-spline  collocation  and  compact  finite 
difference  schemes,  the  two  best  methods  discussed  above. 
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4.  WAVE  EQUATION  IN  BOUNDED  DOMAINS 

In  this  section,  B-spline  collocation  methods  and  compact  finite  difference  methods  are  used 
to  solve  the  wave  equation  in  non-periodic  domains.  The  problem  is  defined  as 


ut  +  ux  =  0  for  0  <  x  <  1, 

u(0,  t)  —  exp  (-ikt)  (10) 

The  exact  solution  assuming  periodicity  in  time  is 

w(x,  t)  =  exp  (ik(x  —  t))  (11) 

For  numerical  solution,  it  is  assumed  that  u(x,t)  takes  the  form  u(x,t)  =  v(x)exp(—ikt)  and 
solve  the  following  equations  for  v(x)  : 

dv 

ikv  —  — , 
ax 

v(0)  =  1  (12) 

The  equation  is  discretized  with  B-spline  collocation  and  compact  finite  difference  schemes  on 
both  uniform  and  non-uniform  grids. 


4.1.  B-spline  Collocation  Formulation 

As  mentioned  in  section  2.1.3,  collocation  points  at  the  B-spline  maxima  are  selected.  In 
general,  using  this  “B-spline  maxima”  collocation  formulation  with  splines  of  order  d,  matrices 
with  d  - 1-1  non-zero  diagonals  will  be  obtained.  In  the  case  of  uniform  grids  away  from  the 
boundary,  there  are  only  d  non-zero  diagonals  as  the  maxima  of  splines  coincide  with  the  knot 
points.  After  discretization,  a  matrix  equation  iuMa  =  Dia  is  obtained,  where  M  and  D\  are 
the  mass  and  first  derivative  operator  matrix,  and  a  is  the  B-spline  coefficient  vector. 

The  boundary  condition  is  implemented  by  replacing  the  operator  at  the  boundary  collocation 
point  with  Vo  =  1. 


4.2.  Compact  Finite  Difference  Formulation 

Lele  presents  a  comprehensive  study  of  high  resolution  finite  difference  schemes  [23].  In  his 
paper,  the  effective  wavenumber  in  a  periodic  domain  is  investigated.  For  domains  with  non¬ 
periodic  boundaries,  the  same  analysis  is  used  to  obtain  the  effective  wavenumbers  both  for  the 
interior  and  the  special  boundary  schemes.  The  effective  wavenumbers  for  the  boundary  schemes 
are  in  general  complex,  with  the  real  part  associated  with  the  dispersive  error  and  the  imaginary 
part  associated  with  the  dissipative  error.  The  conservative  formulation,  eigenvalue  analysis  and 
stability  limits  for  explicit  schemes  are  also  presented.  For  the  details,  the  reader  is  directed 
to  [23]. 

In  this  section,  two  issues  are  addressed.  The  first  is  an  alternative  approach  to  studying  the 
boundary  formulation,  instead  of  the  effective  wavenumber  analysis  of  Lele.  The  second  is  the 
formulation  of  schemes  with  non-uniform  grids.  The  same  problem  is  then  solved  which  offers 
direct  comparison  with  the  B-spline  collocation  method. 

4-2.1.  Boundary  formulation 

To  discretize  the  hyperbolic  equation,  the  numerical  schemes  need  to  resolve  the  traveling  waves 
in  the  domain.  The  boundary  formulation  is  studied  using  normal  modal  analysis.  Normal  modal 
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analysis  is  also  used  by  Carpenter  et  al  [5]  to  investigate  the  stability  of  boundary  treatments 
for  compact  finite  difference  schemes.  The  similarities  of  these  two  analyses  will  be  pointed  out 
after  the  description  of  the  current  boundary  formulation. 

In  the  interior,  the  compact  finite  difference  approximation  of  the  first  derivative  is.  derived 
from  equation  (8),  which  can  be  rewritten  more  generally  as  : 


v'i+Y,ai(vi+j+v'i-j)  =  Y,a- 

j= i  i 


Vj+j  Vj—j 
'*  2jA.x 


(13) 


where  m  is  related  to  the  matrix  bandwidth  w  by  m  =  3£rp'.  Knowing  that  v'  =  ikv,  the 
following  expression  is  obtained  : 

m  m 

X]  civi-j  +  ikVi  -  XZ  cj*vi+j  =  0  (14) 

3=1  3=1 


where  cj  =  +ictjkAx  and  Cj*  is  the  conjugate  of  Cj.  This  can  be  interpreted  as  a  linear  recur¬ 
sion  relation  for  u*.  Such  a  recursion  has  solutions  AJ,  where  A  is  a  function  of  fc.  Substituting 
Vj  =  AJ  into  equation  (14),  the  characteristic  polynomial  is  obtained  : 

m  m 

X>A-’+ifc-X>*Aj=0  (15) 

j=l  j= 1 

which  has  2 m  roots.  In  general,  if  A+  is  a  root,  then  A_  =  A+*-1  is  also  a  root.  These  root 
pairs  are  denoted  as  type  I  root  pairs.  If  |A|  =  1,  A  =  A*-1.  In  this  case  there  can  be  two 
independent  roots.  These  roots  are  denoted  as  type  II  roots.  In  the  limit  k  ->  0,  equation  (15) 
has  the  form 


tt* 

£g(A-J-A')  =  0 


i= 1 


(16) 


which  always  has  the  solutions 


A  -  A-1  =  0  A  =  ±1 


(17) 


Changing  notation  to  that  of  effective  wavenumbers, 


A  =  exp(ikAx)  =>•  u(x ,  t)  =  exp 


(18) 


type  I  root  pairs  correspond  to  conjugate  pairs  of  complex  effective  wavenumbers  k  and  k* ,  while 
type  II  roots  yield  real  k.  Conjugate  pairs  of  complex  effective  wavenumbers  represent  a  pair 
of  solutions,  one  of  which  grows  exponentially  in  amplitude  to  the  right,  the  other  to  the  left. 
Also,  for  k  =  0,  the  two  solutions  yield  k  =  0  and  k  =  kmax.  Clearly,  of  the  2 m  solutions, 
only  one  solution  with  real  k  can  be  a  valid  approximation  to  the  exact  solution.  The  remainder 
are  spurious.  When  equation  (13)  is  used  to  solve  equation  (12),  the  coefficients  of  the  various 
solutions  are  determined  by  the  boundary  conditions  and  special  differencing  schemes  used  near 
the  boundaries.  Clearly,  the  boundary  schemes  should  be  chosen  to  make  the  amplitudes  of 
spurious  solutions  as  small  as  possible. 
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To  see  how  this  works,  consider  the  tridiagonal  and  pentadiagonal  interior  scheme  (see  table 
1).  For  these  two  cases,  the  coefficients  in  the  characteristic  polynomials  and  their  corresponding 
roots  are  given  in  table  5.  A2  and  A3  are  complex  conjugate  pairs  while  Ao  and  Ax  have  magnitude 
1.  A0  represents  the  approximation  to  the  exact  solution  exp(ikAx)  to  the  order  associated  with 
the  scheme  and  it  has  a  positive  group  velocity.  Ai  is  a  spurious  wave  with  a  negative  group 
velocity.  A2  and  A3  are  spurious  waves  growing  exponentially  in  magnitude  to  the  right  and  left 
respectively.  For  the  spurious  waves  that  grow  exponentially  to  the  right,  the  magnitude  of  the 
waves  is  largest  at  the  right  boundary.  Thus,  by  arranging  the  right  boundary  schemes  to  make 
the  A 2  wave  (for  example)  small  at  the  right  boundary,  the  A2  solution  is  small  everywhere, 
regardless  of  the  length  of  the  domain.  Similarly  waves  growing  to  the  left  (e.g.  A3),  should 
be  controlled  at  the  left  boundary.  For  waves  with  |A|  =  1,  or  equivalently  real  fc,  the  “group 
velocity”  vg  =  dk/dk  determines  which  boundary  should  control  the  wave.  With  positive  group 
velocity,  the  left  boundary  controls  the  wave  because  when  solving  the  transient  problem  (10), 
information  from  the  boundary  will  propagate  into  the  domain  from  the  left.  Similarly,  waves 
with  negative  group  velocity  are  controlled  at  the  right  boundary.  Thus  the  spurious  solution  Ai 
will  be  controlled  by  the  right  boundary  scheme,  while  the  physical  boundary  conditions  at  the 
left  boundary  will  control  the  physical  solution  Aq. 

To  determine  the  appropriate  inflow  boundary  schemes,  consider  the  general  solution,  which 
near  the  inflow  boundary  can  be  written  as 

Vj  =  p0Aj)  +  p3A33  +  0((kAx)n)  (19) 

where  n  is  the  order  of  the  error  in  the  interior  scheme  (5  or  9  for  tridiagonal  and  pentadiagonal 
schemes  respectively).  Note  that  for  the  tridiagonal  scheme,  A2  and  A3  can  be  considered  to  be 
zero.  The  0((kAx)n)  term  is  the  contribution  of  the  Ai  and  A2  waves,  which  will  be  this  small  by 
construction  of  the  right  boundary  schemes.  Using  this  expression,  the  left  boundary  schemes  are 
constructed  to  make  p0  =  1  +  0((fcA:r)n)  and  ps  =  0((kAx)n)  (for  the  pentadiagonal  scheme). 
This  is  accomplished  using  schemes  of  the  form 


m+i  ^  3m— i 

y  ctijVj  =  &ijVj  for  0  <  i  <  m  —  1 

j= 0  j= 0 


(20) 


for  the  first  m  =  points,  except  for  the  boundary  point  (i  =  0)  which  is  replaced  by  the 
boundary  condition  Vo  =  1.  The  coefficients  for  bandwidth  3  and  5  (m  =  1  and  2)  are  shown 
in  table  6.  A  Taylor  series  analysis  of  these  schemes  shows  them  to  be  of  the  same  order  as 
the  interior  scheme,  and  indeed  this  is  how  they  were  derived.  This  appears  to  be  a  sufficient 
condition  for  the  suppression  of  the  spurious  waves  to  the  desired  order,  though  the  theory  of 
Gustafsson  ([15]  and  [5])  implies  that  boundary  schemes  one  order  lower  that  the  interior  should 
be  adequate. 

Near  the  outflow  boundary,  the  general  solution  can  be  written  similarly  to  (19) 

VN-I  =  P0A0'  +  pi  Af '  +  P2A21  +  0({kAx)n)  (21) 


where  p\  =  p\A^  and  N  is  the  grid  number  of  the  right  boundary.  Boundary  schemes  that  are 
“mirror  images”  of  the  left  boundary  scheme  result  in  p[  and  pf2  =  0((fcA:r)n)  (pentadiagonal 
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TABLE  5 

Coefficients  and  roots  of  characteristic  polynomials.  The  various  coefficients 
in  the  expressions  for  A  are  given  to  four  digit  accuracy. 

TABLE  5a 

Coefficients  of  characteristic  polynomials 


Band¬ 

width 

Cl 

C2 

3 

5 

|  +  \ikAx 
§  +  %ikAx 

N.A. 

+  4zikAx 

216  *  36 

TABLE  5b 

Roots  of  characteristic  polynomials 

Band¬ 

width 

Ao 

Aj 

A2 

A3 

3 

exp(ifcAx) 

+0((*A*)5) 

-1.0000 

+0.3333  ikAx 
+0.0555(fcAx)2  +  •  •  * 

N.A. 

N.A. 

5 

exp(ifcAa:) 

+0((kAxf) 

-1.0000 

+0.1636  ikAx 
+0.0134(fcAx)2  +  •  ■  • 

-6.2397 

+1.1118  ikAx 
-0.0733(fcAz)2  +  ■  ■  ■ 

-0.1603 

+0.0286  ikAx 
+0.0070(fcAz)2  +  •  •  • 

scheme).  Thus  we  have 


m+N—i  j  3m— N+i 

y,  OLN-i  jv'N_j  =  —  y  -a,N-i  jVN-j  for  N>i>N  —  m+l  (22) 
j= o  j— 0 

where  again  the  coefficients  are  given  in  table  6. 

The  boundary  scheme  analysis  presented  here  is  similar  to  the  GKS  stability  analysis  of  bound¬ 
ary  treatments  in  Carpenter  et  al  [5],  in  which  a  similar  model  problem  is  used  and  in  which 
the  same  spurious  waves  are  treated.  However,  in  Carpenter  et  al ,  the  assumed  temporal  form  of 
the  solution  is  more  general  in  that  the  frequency  k  (in  our  notation,  see  (10))  is  allowed  to  be 
complex.  The  concern  is  then  whether  the  time-periodic  solutions  of  (10)  of  the  form  used  here 
are  stable.  For  the  fourth  order  tridiagonal  scheme,  Carpenter  et  al  show  the  combined  interior 
and  boundary  schemes  to  be  GKS  stable,  but  they  do  not  treat  the  eighth  order  pentadiagonal 
scheme  discussed  here.  It  seems  plausible  that  interior  and  boundary  schemes  that  consistently 
suppress  the  spurious  waves  as  discussed  here  will  be  GKS  stable.  But  this  remains  to  be  seen. 


4-2.2.  Non-uniform  grids 

Another  issue  that  needs  to  be  addressed  is  the  formulation  of  the  compact  finite  difference 
scheme  with  non-uniform  grids.  The  approach  is  to  apply  a  mapping  which  uses  the  uniform  mesh 
scheme  in  the  mapped  coordinate.  The  mesh  mapping  is  given  in  equation  (7).  The  discretization 
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TABLE  6 


Coefficients  of  the  boundary  formulation  for  the  first  derivative. 


Band- 

i 

on 

OLi2 

a, 3 

width 

3 

0 

1 

3 

0 

0 

5 

0 

1 

12 

15 

0 

5 

1 

1 

15 

1 

2 

2 

3 

Band¬ 

width 

i 

OtO 

an 

fl<2 

<2i3 

di4 

a*5 

&i6 

3 

0 

17 

6 

3 

2 

3 

2 

6 

0 

0 

0 

K 

n 

79 

77 

55 

20 

5 

1 

1 

0 

u 

20 

5 

4 

3 

4 

5 

60 

r; 

1 

247 

19 

1 

13 

1 

1 

n 

0 

1 

900 

12 

3 

9 

12 

300 

u 

equations  (8)  ,(20)  and  (22)  are  modified  by  the  mesh  mapping 


dv 

dv 

di 

dx 

=  di 

i  dx 

Thus  the  equation  to  be  solved  is  =  ikv .  The  same  interior  and  boundary  scheme  are  then 
used  in  f . 

4.3.  Comparison 

Tests  based  on  the  solution  of  the  wave  equation  were  carried  out  with  N  =  100.  Before 
discussing  the  results,  however,  it  should  be  noted  that  different  from  the  effective  wavenumbers, 
the  accuracy  of  the  solution  of  the  wave  equation  is  dependent  of  the  number  of  intervals  N  apart 
from  the  wavenumber.  In  this  sense  the  results  here  are  less  general  than  those  of  the  effective 
wavenumber.  Nevertheless,  using  the  same  N  for  both  schemes  allows  us  to  compare  their  order 
of  convergence  and  resolution. 

The  results  on  uniform  grids  are  discussed  first.  The  L 2  errors  in  the  representation  of  the 
solution  of  the  wave  problem  using  B-spline  collocation  and  compact  finite  difference  methods  are 
shown  in  figure  7.  For  B-spline  collocation  methods,  the  errors  vary  with  k  like  kw+2,  Notice  that 
in  periodic  domains,  the  convergence  rates  of  the  eigenvalue  of  the  first  derivative  operator  and 
the  eigenfunction  are  kw+2  and  kw+l  respectively  (see  figures  4,  6  and  table  3).  For  compact  finite 
difference  schemes,  the  L2  error  varies  with  k  like  k2w~l .  This  is  consistent  with  the  theoretical 
convergence  rate  (figure  4  and  table  3).  However,  the  curve  is  not  smooth  as  there  are  many 
small  wriggles  on  it.  Apparently,  the  boundary  condition  and  boundary  schemes  do  not  affect 
their  convergence  rate  of  either  scheme. 

As  mentioned  before,  changing  N  would  affect  the  error  curves.  In  particular,  the  error  curves 
shift  vertically.  This  is  because  for  both  schemes,  the  local  truncation  error  and  hence  the  L2 
error  have  the  form  a(Ax)l~l  f^(x),  where  a  is  a  real  constant  and  i  is  an  integer  ([23]  and  [31]). 
One  way  to  obtain  a  curve  that  is  valid  for  all  N  is  to  plot  error /N  versus  k/kmax.  The  resulting 
curve  would  not  shift  as  N  changes  except  when  the  error  is  close  to  1. 
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fc/^max  fc/fcmax 


FIG.  7.  Z/2  error  in  the  representation  of  the  solution  of  the  wave  equation  with  wavenumber  k  on  uniform 

grids  for  matrix  bandwidth  (a) 3,  (b) 5:  - ,  b-spline  collocation; - ,  compact  finite  difference. 


Perhaps  more  important  than  the  order  of  convergence  is  the  resolution  of  the  two  schemes. 
The  well-resolved  fraction  of  the  wavenumber  range  for  the  solution  of  the  wave  equation  is  shown 
in  table  7.  It  can  be  seen  that  for  tridiagonal  schemes,  the  two  have  almost  the  same  resolution. 
For  pentadiagonal  schemes,  compact  finite  difference  has  better  resolution  due  to  the  higher  order 
of  convergence. 

Another  issue  is  the  effect  of  a  non-uniform  grid.  The  L2  errors  in  the  representation  of 
the  solution  of  the  wave  problem  using  the  two  numerical  schemes  on  non-uniform  grids  are 
shown  in  figure  8.  Here,  the  wavenumber  is  normalized  by  the  maximum  wavenumber  kmax  = 
-g— — ,  where  Axmax  is  the  maximum  grid  spacing.  Basically,  both  schemes  maintain  the  same 
convergence  rate  as  in  the  case  of  uniform  grids.  Note  that  for  the  compact  finite  difference 
schemes,  the  curves  turn  up  at  the  lowest  wave  number  and  the  cause  is  not  clear.  With  regard 
to  the  resolved  fraction,  table  8  indicates  that  the  two  tridiagonal  schemes  again  have  the  same 
resolution.  (Note  however  that  in  non-uniform  grids,  B-spline  collocation  has  elements  outside 
the  three  “main”  diagonals  in  the  interior.)  For  pentadiagonal  schemes,  compact  finite  difference 
has  better  resolution. 


The  order  of  convergence  of  the  two  schemes  suggests  that  the  difference  in  resolution  properties 
between  compact  finite  difference  and  B-spline  collocation  will  become  bigger  as  the  matrix 
bandwidth  increases.  However,  the  wriggles  on  the  compact  finite  difference  curves  may  decrease 
its  resolution.  Also,  comparing  results  in  periodic  domains  (table  7  and  table  4),  it  is  found 
that  the  resolution  in  finite  domains  is  substantially  lower.  In  particular,  the  error  reaches  1  at 
k/kmax  &  0.6  and  then  levels  off. 

Comparison  of  the  numerical  solution  and  the  analytical  solution  as  functions  of  x  suggests 
what  happens  at  the  plateau  of  the  error  curve.  Before  the  error  plateau,  the  numerical  solution  is 
capturing  the  propagating  wave,  although  the  numerical  and  analytical  solution  are  getting  more 
and  more  out  of  phase  as  x  increases  in  the  domain.  This  is  expected  as  the  effective  wavenumber 
is  not  exact.  At  the  error  plateau,  the  numerical  solution  cannot  capture  the  propagating  wave 
at  all.  Small  wriggles  appear  at  the  inflow  boundary  and  gradually  die  away  as  x  increases. 
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FIG.  8.  L2  error  in  the  representation  of  the  solution  of  the  wave  equation  with  wavenumber  k  on  non- 
uniform  grids  for  matrix  bandwidth  (a)3,  (6)5:  - ,  b-spline  collocation; - ,  compact  finite  difference. 


TABLE  7 

Resolved  fraction  for  the  solution  of  the  wave  equation  for  uniform  grid  distribution 


Bandwidth  B-spline  collocation  compact  finite  difference 


10  % 

1  % 

0.1  % 

10  % 

1  % 

0.1  % 

3 

0.25 

0.15 

0.10 

0.24 

0.15 

0.09 

5 

0.40 

0.30 

0.22 

0.45 

0.34 

0.27 

TABLE  8 

Resolved  fraction  for  the  solution  of  the  wave  equation  for  non-uniform  grid  distribution 


Bandwidth  B-spline  collocation  compact  finite  difference 


10  % 

1  % 

0.1  % 

10  % 

1  % 

0.1  % 

3 

0.29 

0.19 

0.12 

0.29 

mm 

0.12 

5 

0.47 

0.35 

0.25 

EEB 

0.34 
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5.  HEAT  EQUATION  IN  BOUNDED  DOMAINS 

In  this  section,  the  eigenvalue  problem  arising  from  the  heat  equation  is  solved  using  B-spline 
collocation  and  compact  finite  difference  methods.  The  problem  is  defined  as 


vn  =  Xv  for  0  <  x  <  1,  (24) 

with  some  boundary  conditions,  the  most  common  ones  being  the  Dirichlet  (ir(0)  =  v(l)  —  0) 
and  Neumann  (t/'(0)  =  z;'(l)  =  0)  boundary  conditions.  In  both  cases,  the  eigenvalues  are 

\k  =  -{irk)2,  (25) 

where  k  is  an  integer.  The  corresponding  eigenfunctions  are  Vk{x)  =  sin(7r kx)  and  ^(x)  = 
cos(7rfcx)  for  Dirichlet  and  Neumann  boundary  conditions  respectively. 

5.1.  B-spline  Collocation  Formulation 

Discretizing  with  B-spline  collocation  method,  we  obtain  the  matrix  equation  A Ma  =  D2a, 
where  M  is  the  mass  and  D2  the  second  derivative  operator  matrix,  and  a  the  B-spline  coefficient 
vector  for  the  eigenfunctions.  For  Dirichlet  boundary  conditions,  v0  =  vn  ^  0.  For  Neumann 
boundary  conditions,  v*0  and  v'N  are  set  to  zero. 

5.2.  Compact  Finite  Difference  Formulation 

Similar  to  the  case  of  first  derivative,  the  discretized  derivative  operators  are  derived  from 
equation  (9)  in  the  interior.  Near  the  boundary,  the  symmetry  breaks  down  and  the  corresponding 
equation  becomes 


m+i  ..  3m+l- i 

ai3vj  =  13  aiivi  for  o  <  *  <  m  —  1,  (26) 

j= 0  ^  '  j= 0 

m+N~i  ^  3m+l— JV+i 

aN-i  jv'n-j  =  >2  Y!  iVN~3  iorN>i>N  —  m  +  l 

j=0  ^  X'  j- 0 

where  m  =  The  coefficients  in  equation  (26)  are  determined  by  matching  the  Taylor  series 
coefficients  to  one  order  less  than  the  interior  for  tridiagonal  schemes  and  to  the  same  order  of 
the  interior  for  pentadiagonal  schemes  (using  boundary  of  the  same  order  for  tridiagonal  schemes 
gives  rise  to  bad  resolution  for  unknown  reasons).  The  coefficients  are  shown  in  table  9  for  the 
two  schemes.  After  discretization,  a  generalized  eigenvalue  problem  Ail^a  =  D2ol  is  obtained, 
where  M2  is  the  mass  and  D2  the  second  derivative  operator  matrix,  and  a  the  eigenfunction. 
The  generalized  eigenvalue  problem  can  be  solved  with  the  appropriate  boundary  condition. 

5.2.1.  Boundary  conditions 

The  Dirichlet  boundary  conditions  are  implemented  by  setting  vq  =  vn  =  0.  For  Neumann 
boundary  conditions,  a  one-sided  explicit  (i.e.  not  compact)  finite  difference  scheme  is  used  to  set 
v'0  =  v'N  =  0.  Note  that  this  makes  the  boundary  scheme  inconsistent  with  the  interior  scheme. 
Also,  a  very  long  one-sided  finite  difference  expression  is  required  to  maintain  the  same  order  as 
the  interior  compact  finite  difference  approximations. 
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TABLE  9 

Coefficients  of  the  boundary  formulation  for  the  second  derivative 


3.34302  x  10" 


3.36092  X  101 


1.17128  X  10^ 


8.58850  X  10" 


2.15798  x  101 
6.76347  x  10” 


1.03882  X  10z 


-2.92763  X  10z 


1.91514  x  10z 


-2.87907  X  101 
5.39365  X  10”1 
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5.2.2.  Non-uniform  grids 

To  solve  the  problem  on  a  non-uniform  grid,  a  mesh  mapping  is  used  as  in  the  wave  equation. 
Notice  from  the  chain  rule 

—  =  (  *i\ 2  cP£  dv 

dx2  \dx)  d?  dx2  d£,  1  ' 

The  derivative  in  the  non-uniform  x  coordinate  is  expressed  in  terms  of  those  in  the  transformed 
uniform  £  coordinate.  In  the  £  coordinate,  there  are  finite  difference  representation  of  the  deriva¬ 
tive  operators  (expressed  as  M±lD\  and  Note  that  Mi  and  M2  are  different).  The 

finite  difference  approximation  of  the  second  derivative  operator  can  hence  be  expressed  as  in 
equation  (27). 


5.3.  Comparison 

Tests  based  on  the  eigenvalue  problem  are  performed  using  N  =  30.  Results  based  on  dif¬ 
ferent  N  suggest  that  N  has  no  influence  on  the  order  of  convergence  and  minor  influence  on 
the  well-resolved  fraction.  The  results  obtained  on  uniform  grids  are  presented  first.  The  errors 
in  approximating  the  eigenvalues  are  shown  in  figure  9.  Regardless  of  the  boundary  conditions, 
B-spline  collocation  schemes  have  eigenvalue  errors  which  decrease  with  wavenumber  as  A;™*1, 
while  compact  finite  difference  has  a  convergence  rate  of  k2w.  Both  of  the  above  are  consistent 
with  their  corresponding  convergence  rates  in  periodic  domains  (see  table  3).  For  compact  finite 
difference,  however,  the  boundary  conditions  do  have  an  effect  on  the  magnitude  of  the  error. 
Neumann  boundary  conditions  give  larger  errors  in  the  eigenvalues,  perhaps  due  to  the  boundary 
approximation  of  vf.  Also  with  the  compact  finite  difference,  there  are  some  sharp  decrease  in 
error  at  particular  wavenumbers  for  reasons  that  are  not  clear.  At  high  wavenumbers,  wriggles 
appear  on  the  compact  finite  difference  curves  irrespective  of  the  boundary  conditions.  With 
regard  to  resolution,  we  refer  to  table  10,  which  gives  the  resolved  fraction  for  the  eigenvalues. 
In  many  cases,  compact  finite  difference  schemes  provide  better  resolution  for  the  eigenvalues. 
However,  for  pentadiagonal  scheme,  B-splines  have  better  resolved  fractions  in  many  cases.  How¬ 
ever,  due  to  the  high  convergence  order  of  the  compact  finite  difference,  the  more  stringent  the 
tolerance  for  resolved  fractions,  the  better  the  compact  finite  difference  does. 

The  I/2  errors  of  the  eigenfunctions  of  the  heat  equation  are  shown  in  figure  10.  For  the 
B-spline  collocation  schemes,  the  convergence  rate  for  both  Dirichlet  and  Neumann  boundary 
conditions  appears  to  be  fcw+1,  but  in  the  Neumann  case  this  asymptotic  rate  is  not  attained 
until  k  <  0.06,  with  the  resulting  impact  on  resolution.  For  compact  finite  difference  schemes,  the 
eigenfunctions  have  errors  that  converge  at  a  rate  approximately  equal  to  k2w.  However,  with 
Neumann  boundary  conditions,  the  errors  are  again  larger.  In  fact,  using  Dirichlet  boundary 
condition,  the  tridiagonal  schemes  give  a  solution  that  is  exact  to  round-off  errors.  It  is  also  very 
interesting  to  note  that  the  pentadiagonal  scheme  curve  shows  two  sharp  decrease  at  yf-X /kmax  = 
|  and  y/^X/kmax  —  §.  At  these  two  particular  wavenumbers,  the  symmetries  of  the  approximate 
eigenfunctions  make  the  point  representations  exact.  Table  11  shows  the  resolved  fraction  of 
eigenfunctions.  Same  as  for  the  eigenvalues,  compact  finite  difference  schemes  in  general  provide 
better  resolution  for  tridiagonal  methods  while  B-splines  do  better  for  pentadiagonal  schemes. 

On  non-uniform  grids,  the  behavior  of  both  B-spline  collocation  and  compact  finite  difference 
schemes  is  shown  in  figures  11,  12, tables  12  and  13.  The  errors  in  the  eigenvalues  of  the  heat 
equation  are  shown  in  figure  11.  B-spline  collocation  methods  maintain  the  same  convergence 
rate  of  as  in  the  case  of  uniform  grids  irrespective  of  the  boundary  conditions.  Compact 
finite  difference  schemes,  however,  show  a  degradation.  Convergence  rates  of  the  eigenvalues  is 
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FIG.  9.  Error  of  the  eigenvalues  of  the  heat  equation  with  wavenumber  k  on  uniform  grids  for  matrix 

bandwidth  (a)3,  (6)5:  - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both  for 

b-spline; - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both  for  compact  finite 

difference 


FIG.  10.  Error  of  the  eigenfunctions  of  the  heat  equation  with  wavenumber  k  on  uniform  grids  for 

matrix  bandwidth  (a)3,  (6)5:  - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both 

for  b-spline; - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both  for  compact 

finite  difference.  For  the  compact  finite  difference,  with  bandwidth  of  3  and  Dirichlet  boundary  conditions,  the 
eigenfunctions  are  exact  to  round-off  errors. 
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TABLE  10 

Resolved  fraction  for  the  eigenvalues  for  uniform  grid  distribution 


Bandwidth  Boundary  B-spline  collocation  compact  finite 

condition  difference 


10  % 

i  % 

0.1  % 

10  % 

i  % 

0.1  % 

3 

Dirichlet 

0.33 

0.10 

0.03 

0.66 

0.36 

0.20 

3 

Neumann 

0.36 

0.10 

0.03 

0.46 

0.36 

0.16 

5 

Dirichlet 

0.80 

0.46 

0.26 

>  0.76 

0.50 

0.40 

5 

Neumann 

0.76 

0.46 

0.26 

0.43 

0.36 

0.23 

TABLE  11 

Resolved  fraction  for  the  eigenfunctions  for  uniform  grid  distribution 


Bandwidth  Boundary  B-spline  collocation  compact  finite 
condition  difference 


10  % 

i  % 

0.1  % 

10  % 

i  % 

0.1  % 

3 

Dirichlet 

0.43 

0.23 

0.13 

1.00 

1.00 

1.00 

3 

Neumann 

0.40 

0.13 

0.06 

0.23 

0.16 

0.10 

5 

Dirichlet 

0.73 

0.46 

0.30 

0.43 

0.37 

0.33 

5 

Neumann 

0.93 

0.56 

0.33 

0.36 

0.20 

0.16 

26 


W.Y.Kwok,  R.D. Moser,  and  J. Jimenez 


FIG.  11.  Error  of  the  eigenvalues  of  the  heat  equation  with  wavenumber  k  on  non-uniform  grids  for  matrix 

bandwidth  (a)3,  (6)5:  - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both  for 

b-spline; - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both  for  compact  finite 

difference 


TABLE  12 

Resolved  fraction  for  the  eigenvalues  for  non-uniform  grid  distribution 


Bandwidth  Boundary  B-spline  collocation  compact  finite 
condition  difference 


10% 

i  % 

0.1  % 

10% 

i  % 

0.1  % 

3 

Dirichlet 

0.44 

0.09 

<  0.05 

0.73 

0.44 

0.19 

3 

Neumann 

0.44 

0.14 

0.04 

0.74 

0.54 

0.44 

5 

Dirichlet 

0.92 

0.53 

0.29 

0.73 

0.49 

0.34 

5 

Neumann 

0.97 

0.53 

0.29 

0.54 

0.29 

0.14 

2  to  3  orders  less  than  the  corresponding  rate  of  k2w  on  uniform  grids,  with  Neumann  boundary 
conditions  giving  worse  convergence  rates.  Regarding  resolution,  compact  finite  difference  pro¬ 
vides  better  resolution  for  bandwidth  w  =  3  while  B-spline  collocation  schemes  provides  better 
resolution  for  bandwidth  w  =  5. 

The  Z/2  errors  of  the  eigenfunctions  of  the  heat  equation  are  shown  in  figure  12.  B-spline 
collocation  schemes  give  convergence  rates  of  kw  approximately,  with  Dirichlet  boundary  condi¬ 
tions  giving  slightly  better  solutions  at  low  k.  The  degradation  of  resolution  is  not  serious  when 
non-uniform  grids  are  used  instead  of  uniform  ones.  Compact  finite  difference  schemes,  however, 
show  quite  serious  degradation  of  convergence  and  resolution  on  non-uniform  grids.  They  have 
convergence  rates  of  about  kw,  compared  to  k2w  on  uniform  grids.  A  very  interesting  result  is 
that  B-spline  and  compact  finite  difference  schemes  appear  to  have  the  same  convergence  rates 
on  non-uniform  grids.  From  table  13,  it  can  also  be  seen  that  B-spline  collocation  methods  have 
better  resolution  properties  on  non-uniform  grids. 
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FIG.  12.  Error  of  the  eigenfunctions  of  the  heat  equation  with  wavenumber  k  on  non-uniform  grids  for 

matrix  bandwidth  (a) 3,  (6)5:  - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both 

for  b-spline; - ,  Dirichlet  boundary  condition; - ,  Neumann  boundary  condition,  both  for  compact  finite 

difference. 


TABLE  13 


Resolved  fraction  for  the  eigenfunctions  for  non-uniform  grid  distribution 


Bandwidth  Boundary  B-spline  collocation  compact  finite 
condition  difference 


10  % 

i  % 

0.1  % 

10  % 

i  % 

0.1  % 

3 

Dirichlet 

0.39 

0.19 

0.09 

0.44 

0.24 

0.14 

3 

Neumann 

0.39 

0.14 

0.04 

0.29 

0.09 

<0.04 

5 

Dirichlet 

0.72 

0.43 

0.29 

0.39 

0.29 

0.24 

5 

Neumann 

0.72 

0.43 

0.29 

0.25 

0.15 

0.09 
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6.  DISCUSSION  AND  CONCLUSIONS 

The  results  of  this  paper  indicate  that  in  many  situations  compact  finite  difference  schemes 
have  better  resolution  and  convergence  properties  than  the  other  numerical  methods  tested. 
The  comparisons  were  done  for  schemes  with  the  same  matrix  bandwidths,  which  we  use  as  a 
surrogate  for  computational  cost.  Furthermore,  it  was  shown  that  finite  element  and  B-spline 
Galerkin  methods  had  inferior  resolution  to  compact  finite  difference  and  B-spline  collocation. 
There  are  several  aspects  of  these  results  that  deserve  further  discussion. 

Regarding  high-order  finite  element  methods,  it  was  noted  that  a  reason  for  their  relatively 
poor  performance  in  these  tests  was  their  low-order  (Co)  continuity  at  the  element  boundaries 
(i.e.  knots),  whereas  the  spline  basis  retains  as  high  a  degree  of  continuity  as  possible,  given 
the  order  of  the  piecewise  polynomial  representation.  In  essence,  in  spline  methods,  an  increase 
in  the  degree  of  the  polynomials  is  used  to  increase  the  degree  of  continuity,  while  in  Co  finite 
elements,  it  is  used  to  increase  the  number  of  degrees  of  freedom  of  the  representation.  The 
results  of  the  tests  here  suggest  that  the  added  degrees  of  freedom  do  not  produce  much  in  the 
way  of  added  accurately  represented  modes,  resulting  in  poor  resolution  properties.  However,  the 
improved  resolution  of  splines  is  not  without  cost;  that  is,  the  representation  of  the  polynomials 
in  each  interval  (element)  is  not  isoparametric,  a  very  convenient  property  of  finite  element 
representations.  Consequently,  it  is  much  easier  to  formulate  multi-dimensional  finite  elements 
on  complex  and/or  unstructured  grids,  than  it  is  to  formulate  spline  methods. 

It  was  also  noted  that  piecewise  polynomial  Galerkin  methods  yielded  poorer  representations 
of  complex  exponentials  (the  derivative  eigenfunction)  than  collocation  methods.  This  is  true 
for  both  finite  element  methods  and  spline  methods.  This  is  a  curious  result  because  Galerkin 
approximations  minimize  L<i  error  for  a  given  representation.  The  reason  for  the  curious  result 
is  that  we  are  comparing  methods  with  the  same  matrix  bandwidth.  For  example,  a  Galerkin 
method  that  yields  pentadiagonal  matrices  has  cubic  polynomials,  whereas  a  pentadiagonal  col¬ 
location  methods  has  quintic  polynomials.  The  result  is  a  fourth  order  accurate  representation 
for  Galerkin  and  a  sixth  order  accurate  representation  for  collocation.  However,  there  are  other 
reasons  one  might  choose  a  Galerkin  method,  despite  this  shortcoming;  for  example,  a  Galerkin 
method  is  trivially  shown  to  be  conservative. 

The  two  methods  discussed  here  with  the  best  convergence  and  resolution  properties  are  com¬ 
pact  finite  difference  and  spline  collocation,  and  the  comparison  between  them  includes  4  major 
issues  that  must  be  traded  off  against  the  improved  order  of  accuracy  and  in  many  cases  better 
resolution  of  the  finite  difference  methods: 

1.  The  generally  superior  convergence  and  resolution  of  compact  finite  difference  compared  to 
B-spline  collocation  is  simply  due  to  the  fact  that  in  the  finite  difference  case,  the  “mass  matrix” 
is  not  constrained  to  be  the  same  for  all  derivatives.  There  may,  however,  be  costs  in  code 
complexity  or  computational  effort  in  having  different  mass  matrices,  depending  on  the  details 
of  the  problem  being  solved. 

2.  Another  difference  is  in  the  treatment  near  a  boundary.  In  the  finite  difference  case,  special 
difference  schemes  must  be  formulated  near  the  boundary,  and  such  boundary  treatments  can 
cause  difficulties.  For  the  wave  equation,  a  criterion  for  a  good  boundary  treatment  was  developed 
in  section  4.2.1,  and  schemes  that  satisfy  the  criterion  were  found  by  imposing  a  formal  order 
of  accuracy  consistent  with  the  interior  scheme.  However,  consistent  order  of  accuracy  is  a 
necessary  but  not  necessarily  sufficient  condition  for  the  criterion  to  be  satisfied,  and  directly 
constructing  schemes  to  satisfy  the  criterion  is  prohibitively  cumbersome  in  all  but  the  simplest 
cases  (e.g.  the  tridiagonal  scheme).  Thus  we  do  not  have  a  practical  constructive  prescription 
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for  boundary  schemes  that  satisfy  the  criterion  in  section  4.2.1.  In  the  heat  equation  problem, 
the  development  of  a  similar  criterion  for  good  boundary  schemes  is  not  as  obvious,  so  we  have 
an  even  less  well  defined  procedure  for  the  boundary  treatment  in  this  case.  Finally,  recall  that 
for  the  heat  equation,  with  Neumann  conditions,  an  approximation  of  the  first  derivative  at  the 
boundary  had  to  be  devised,  so  the  derivative  boundary  condition  could  be  imposed.  This  was 
essentially  ad  hoc,  and  was  not  inherently  consistent  with  the  remainder  of  the  scheme.  All  of 
these  difficulties  are  obviated  with  the  spline  method,  since  there  is  no  special  treatment  near 
the  boundary. 

3.  On  a  nonuniform  mesh,  the  spline  method  can  be  used  directly,  without  recourse  to  a 
mapping  to  a  domain  with  a  uniform  mesh,  as  we  did  for  the  finite  difference  case.  Thus  the 
method  can  easily  be  applied  to  an  arbitrary  mesh,  for  which  no  analytic  mapping  is  known. 
Besides,  the  resulting  approximations  are  simpler,  with  no  explicit  metric  terms,  and  in  the  case 
of  approximating  the  second  derivative,  no  first  derivative  term  appears.  Of  course,  one  can 
construct  finite  difference  methods  directly  on  a  nonuniform  mesh  as  well,  but  the  process  is 
cumbersome,  and  generally  yields  schemes  that  are  lower  order  than  the  uniform  mesh  schemes 
(for  the  same  matrix  bandwidth). 

4.  In  sections  4  and  5,  the  error  associated  with  the  spline  collocation  method  was  well-behaved 
and  consistent  with  the  results  of  the  Fourier  analysis  in  Section  3.  The  same  cannot  be  said  for 
the  finite  difference  schemes.  For  them,  the  error  spectra  were  more  erratic,  with  a  variety  of 
unexplained  features.  Furthermore,  in  at  least  one  case  (i.e.  nonuniform  mesh  heat  equation  with 
Neumann  conditions),  the  convergence  rates  appeared  to  be  the  same  as  its  spline  counterpart, 
inconsistent  with  the  simple  Fourier  analysis  in  Section  3. 

Thus,  when  using  a  high-order  spline  collocation  scheme  instead  of  compact  finite  difference 
with  the  same  bandwidth,  one  is  trading  away  a  potentially  higher  convergence  rate  and  somewhat 
better  resolution  in  many  cases  for  a  more  straightforward  and  robust  formulation.  And,  as 
suggested  by  the  results  of  section  5,  in  complicated  situations,  there  is  no  guarantee  that  the 
finite  difference  method  would  actually  yield  the  theoretical  higher  convergence  rate. 
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